/[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

trunk/escript/src/DataLazy.cpp revision 2721 by jfenwick, Fri Oct 16 05:40:12 2009 UTC trunk/escriptcore/src/DataLazy.cpp revision 6168 by jfenwick, Wed Apr 13 03:08:12 2016 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2016 by The University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Apache License, version 2.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.apache.org/licenses/LICENSE-2.0
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14    *
15    *****************************************************************************/
16    
17  #include "DataLazy.h"  #include "DataLazy.h"
 #ifdef USE_NETCDF  
 #include <netcdfcpp.h>  
 #endif  
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
 #ifdef _OPENMP  
 #include <omp.h>  
 #endif  
 #include "FunctionSpace.h"  
 #include "DataTypes.h"  
18  #include "Data.h"  #include "Data.h"
19  #include "UnaryFuncs.h"     // for escript::fsign  #include "DataTypes.h"
20    #include "EscriptParams.h"
21    #include "FunctionSpace.h"
22  #include "Utils.h"  #include "Utils.h"
23    #include "DataVectorOps.h"
24    
25  #include "EscriptParams.h"  #include <iomanip> // for some fancy formatting in debug
26    
27    using namespace escript::DataTypes;
28    
29  #include <iomanip>      // for some fancy formatting in debug  #define NO_ARG
30    
31  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
32  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
# Line 44  bool privdebug=false; Line 40  bool privdebug=false;
40    
41  // #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();}  // #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();}
42    
43  #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}  // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "SIZE LIMIT EXCEEDED " << m_height << endl;resolveToIdentity();}
44    
45    
46    #define SIZELIMIT if (m_height>escript::escriptParams.getTOO_MANY_LEVELS())  {if (escript::escriptParams.getLAZY_VERBOSE()){cerr << "SIZE LIMIT EXCEEDED height=" << m_height << endl;}resolveToIdentity();}
47    
48  /*  /*
49  How does DataLazy work?  How does DataLazy work?
50  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 58  A special operation, IDENTITY, stores an Line 56  A special operation, IDENTITY, stores an
56  This means that all "internal" nodes in the structure are instances of DataLazy.  This means that all "internal" nodes in the structure are instances of DataLazy.
57    
58  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
59  Note that IDENITY is not considered a unary operation.  Note that IDENTITY is not considered a unary operation.
60    
61  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
62  It must however form a DAG (directed acyclic graph).  It must however form a DAG (directed acyclic graph).
# Line 66  I will refer to individual DataLazy obje Line 64  I will refer to individual DataLazy obje
64    
65  Each node also stores:  Each node also stores:
66  - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was  - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
67      evaluated.          evaluated.
68  - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to  - m_buffsrequired ~ the large number of samples which would need to be kept simultaneously in order to
69      evaluate the expression.          evaluate the expression.
70  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
71    
72  When a new node is created, the above values are computed based on the values in the child nodes.  When a new node is created, the above values are computed based on the values in the child nodes.
# Line 109  namespace escript Line 107  namespace escript
107  namespace  namespace
108  {  {
109    
 enum ES_opgroup  
 {  
    G_UNKNOWN,  
    G_IDENTITY,  
    G_BINARY,        // pointwise operations with two arguments  
    G_UNARY,     // pointwise operations with one argument  
    G_UNARY_P,       // pointwise operations with one argument, requiring a parameter  
    G_NP1OUT,        // non-pointwise op with one output  
    G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter  
    G_TENSORPROD,    // general tensor product  
    G_NP1OUT_2P,     // non-pointwise op with one output requiring two params  
    G_REDUCTION      // non-pointwise unary op with a scalar output  
 };  
110    
111    // enabling this will print out when ever the maximum stacksize used by resolve increases
112    // it assumes _OPENMP is also in use
113    //#define LAZY_STACK_PROF
114    
115    
116    
117  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",  #ifndef _OPENMP
118              "sin","cos","tan",    #ifdef LAZY_STACK_PROF
119              "asin","acos","atan","sinh","cosh","tanh","erf",    #undef LAZY_STACK_PROF
120              "asinh","acosh","atanh",    #endif
121              "log10","log","sign","abs","neg","pos","exp","sqrt",  #endif
122              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",  
123              "symmetric","nonsymmetric",  
124              "prod",  #ifdef LAZY_STACK_PROF
125              "transpose", "trace",  std::vector<void*> stackstart(getNumberOfThreads());
126              "swapaxes",  std::vector<void*> stackend(getNumberOfThreads());
127              "minval", "maxval"};  size_t maxstackuse=0;
128  int ES_opcount=43;  #endif
129  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,  
             G_UNARY,G_UNARY,G_UNARY, //10  
             G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17  
             G_UNARY,G_UNARY,G_UNARY,                    // 20  
             G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 28  
             G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35  
             G_NP1OUT,G_NP1OUT,  
             G_TENSORPROD,  
             G_NP1OUT_P, G_NP1OUT_P,  
             G_NP1OUT_2P,  
             G_REDUCTION, G_REDUCTION};  
 inline  
 ES_opgroup  
 getOpgroup(ES_optype op)  
 {  
   return opgroups[op];  
 }  
130    
131  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
132  FunctionSpace  FunctionSpace
133  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
134  {  {
135      // perhaps this should call interpolate and throw or something?          // perhaps this should call interpolate and throw or something?
136      // maybe we need an interpolate node -          // maybe we need an interpolate node -
137      // 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
138      // programming error exception.          // programming error exception.
139    
140    FunctionSpace l=left->getFunctionSpace();    FunctionSpace l=left->getFunctionSpace();
141    FunctionSpace r=right->getFunctionSpace();    FunctionSpace r=right->getFunctionSpace();
142    if (l!=r)    if (l!=r)
143    {    {
144      if (r.probeInterpolation(l))      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
145        if (res==1)
146      {      {
147      return l;          return l;
148      }      }
149      if (l.probeInterpolation(r))      if (res==-1)
150      {      {
151      return r;          return r;
152      }      }
153      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
154    }    }
# Line 187  resultFS(DataAbstract_ptr left, DataAbst Line 160  resultFS(DataAbstract_ptr left, DataAbst
160  DataTypes::ShapeType  DataTypes::ShapeType
161  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
162  {  {
163      if (left->getShape()!=right->getShape())          if (left->getShape()!=right->getShape())
164      {          {
165        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))            if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
166        {            {
167          throw DataException("Shapes not the name - shapes must match for (point)binary operations.");                  throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
168        }            }
169    
170        if (left->getRank()==0)   // we need to allow scalar * anything            if (left->getRank()==0)       // we need to allow scalar * anything
171        {            {
172          return right->getShape();                  return right->getShape();
173        }            }
174        if (right->getRank()==0)            if (right->getRank()==0)
175        {            {
176          return left->getShape();                  return left->getShape();
177        }            }
178        throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");            throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
179      }          }
180      return left->getShape();          return left->getShape();
181  }  }
182    
183  // return the shape for "op left"  // return the shape for "op left"
# Line 212  resultShape(DataAbstract_ptr left, DataA Line 185  resultShape(DataAbstract_ptr left, DataA
185  DataTypes::ShapeType  DataTypes::ShapeType
186  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
187  {  {
188      switch(op)          switch(op)
189      {          {
190          case TRANS:          case TRANS:
191         {            // for the scoping of variables             {                    // for the scoping of variables
192          const DataTypes::ShapeType& s=left->getShape();                  const DataTypes::ShapeType& s=left->getShape();
193          DataTypes::ShapeType sh;                  DataTypes::ShapeType sh;
194          int rank=left->getRank();                  int rank=left->getRank();
195          if (axis_offset<0 || axis_offset>rank)                  if (axis_offset<0 || axis_offset>rank)
196          {                  {
197                 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);              stringstream e;
198              }              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
199              for (int i=0; i<rank; i++)              throw DataException(e.str());
200          {          }
201             int index = (axis_offset+i)%rank;          for (int i=0; i<rank; i++)
202                 sh.push_back(s[index]); // Append to new shape                  {
203              }                     int index = (axis_offset+i)%rank;
204          return sh;             sh.push_back(s[index]); // Append to new shape
205         }          }
206      break;                  return sh;
207      case TRACE:             }
208         {          break;
209          int rank=left->getRank();          case TRACE:
210          if (rank<2)             {
211          {                  int rank=left->getRank();
212             throw DataException("Trace can only be computed for objects with rank 2 or greater.");                  if (rank<2)
213          }                  {
214          if ((axis_offset>rank-2) || (axis_offset<0))                     throw DataException("Trace can only be computed for objects with rank 2 or greater.");
215          {                  }
216             throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");                  if ((axis_offset>rank-2) || (axis_offset<0))
217          }                  {
218          if (rank==2)                     throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
219          {                  }
220             return DataTypes::scalarShape;                  if (rank==2)
221          }                  {
222          else if (rank==3)                     return DataTypes::scalarShape;
223          {                  }
224             DataTypes::ShapeType sh;                  else if (rank==3)
225                 if (axis_offset==0)                  {
226             {                     DataTypes::ShapeType sh;
227                  sh.push_back(left->getShape()[2]);                     if (axis_offset==0)
228                 }                     {
229                 else     // offset==1                          sh.push_back(left->getShape()[2]);
230             {                     }
231              sh.push_back(left->getShape()[0]);                     else         // offset==1
232                 }                     {
233             return sh;                          sh.push_back(left->getShape()[0]);
234          }                     }
235          else if (rank==4)                     return sh;
236          {                  }
237             DataTypes::ShapeType sh;                  else if (rank==4)
238             const DataTypes::ShapeType& s=left->getShape();                  {
239                 if (axis_offset==0)                     DataTypes::ShapeType sh;
240             {                     const DataTypes::ShapeType& s=left->getShape();
241                  sh.push_back(s[2]);                     if (axis_offset==0)
242                  sh.push_back(s[3]);                     {
243                 }                          sh.push_back(s[2]);
244                 else if (axis_offset==1)                          sh.push_back(s[3]);
245             {                     }
246                  sh.push_back(s[0]);                     else if (axis_offset==1)
247                  sh.push_back(s[3]);                     {
248                 }                          sh.push_back(s[0]);
249             else     // offset==2                          sh.push_back(s[3]);
250             {                     }
251              sh.push_back(s[0]);                     else         // offset==2
252              sh.push_back(s[1]);                     {
253             }                          sh.push_back(s[0]);
254             return sh;                          sh.push_back(s[1]);
255          }                     }
256          else        // unknown rank                     return sh;
257          {                  }
258             throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");                  else            // unknown rank
259          }                  {
260         }                     throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
261      break;                  }
262          default:             }
263      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");          break;
264      }          default:
265            throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
266            }
267  }  }
268    
269  DataTypes::ShapeType  DataTypes::ShapeType
# Line 306  SwapShape(DataAbstract_ptr left, const i Line 281  SwapShape(DataAbstract_ptr left, const i
281          throw DataException("Error - Data::swapaxes argument must have at least rank 2.");          throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
282       }       }
283       if (axis0<0 || axis0>rank-1) {       if (axis0<0 || axis0>rank-1) {
284          throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);          stringstream e;
285            e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
286            throw DataException(e.str());
287       }       }
288       if (axis1<0 || axis1>rank-1) {       if (axis1<0 || axis1>rank-1) {
289           throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);          stringstream e;
290            e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
291            throw DataException(e.str());
292       }       }
293       if (axis0 == axis1) {       if (axis0 == axis1) {
294           throw DataException("Error - Data::swapaxes: axis indices must be different.");           throw DataException("Error - Data::swapaxes: axis indices must be different.");
# Line 340  SwapShape(DataAbstract_ptr left, const i Line 319  SwapShape(DataAbstract_ptr left, const i
319  DataTypes::ShapeType  DataTypes::ShapeType
320  GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)  GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
321  {  {
322                
323    // Get rank and shape of inputs    // Get rank and shape of inputs
324    int rank0 = left->getRank();    int rank0 = left->getRank();
325    int rank1 = right->getRank();    int rank1 = right->getRank();
# Line 349  GTPShape(DataAbstract_ptr left, DataAbst Line 328  GTPShape(DataAbstract_ptr left, DataAbst
328    
329    // Prepare for the loops of the product and verify compatibility of shapes    // Prepare for the loops of the product and verify compatibility of shapes
330    int start0=0, start1=0;    int start0=0, start1=0;
331    if (transpose == 0)       {}    if (transpose == 0)           {}
332    else if (transpose == 1)  { start0 = axis_offset; }    else if (transpose == 1)      { start0 = axis_offset; }
333    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)      { start1 = rank1-axis_offset; }
334    else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }    else                          { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
335    
336    if (rank0<axis_offset)    if (rank0<axis_offset)
337    {    {
338      throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");          throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
339    }    }
340    
341    // Adjust the shapes for transpose    // Adjust the shapes for transpose
342    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather    DataTypes::ShapeType tmpShape0(rank0);        // pre-sizing the vectors rather
343    DataTypes::ShapeType tmpShape1(rank1);    // than using push_back    DataTypes::ShapeType tmpShape1(rank1);        // than using push_back
344    for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }    for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
345    for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }    for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
346    
347    // Prepare for the loops of the product    // Prepare for the loops of the product
348    SL=1, SM=1, SR=1;    SL=1, SM=1, SR=1;
349    for (int i=0; i<rank0-axis_offset; i++)   {    for (int i=0; i<rank0-axis_offset; i++)       {
350      SL *= tmpShape0[i];      SL *= tmpShape0[i];
351    }    }
352    for (int i=rank0-axis_offset; i<rank0; i++)   {    for (int i=rank0-axis_offset; i<rank0; i++)   {
353      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
354        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
355      }      }
356      SM *= tmpShape0[i];      SM *= tmpShape0[i];
357    }    }
358    for (int i=axis_offset; i<rank1; i++)     {    for (int i=axis_offset; i<rank1; i++)         {
359      SR *= tmpShape1[i];      SR *= tmpShape1[i];
360    }    }
361    
362    // Define the shape of the output (rank of shape is the sum of the loop ranges below)    // Define the shape of the output (rank of shape is the sum of the loop ranges below)
363    DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      
364    {         // block to limit the scope of out_index    {                     // block to limit the scope of out_index
365       int out_index=0;       int out_index=0;
366       for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z       for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
367       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
# Line 398  GTPShape(DataAbstract_ptr left, DataAbst Line 377  GTPShape(DataAbstract_ptr left, DataAbst
377    return shape2;    return shape2;
378  }  }
379    
380  // determine the number of samples requires to evaluate an expression combining left and right  }       // end anonymous namespace
 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.  
 // The same goes for G_TENSORPROD  
 // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.  
 // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined  
 // multiple times  
 int  
 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  
 {  
    switch(getOpgroup(op))  
    {  
    case G_IDENTITY: return 1;  
    case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);  
    case G_REDUCTION:  
    case G_UNARY:  
    case G_UNARY_P: return max(left->getBuffsRequired(),1);  
    case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);  
    case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);  
    case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);  
    case G_NP1OUT_2P: return 1+max(left->getBuffsRequired(),1);  
    default:  
     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");  
    }  
 }  
   
   
 }   // end anonymous namespace  
   
   
   
 // Return a string representing the operation  
 const std::string&  
 opToString(ES_optype op)  
 {  
   if (op<0 || op>=ES_opcount)  
   {  
     op=UNKNOWNOP;  
   }  
   return ES_opstrings[op];  
 }  
381    
 #ifdef LAZY_NODE_STORAGE  
382  void DataLazy::LazyNodeSetup()  void DataLazy::LazyNodeSetup()
383  {  {
384  #ifdef _OPENMP  #ifdef _OPENMP
# Line 456  void DataLazy::LazyNodeSetup() Line 395  void DataLazy::LazyNodeSetup()
395      m_sampleids[0]=-1;      m_sampleids[0]=-1;
396  #endif  // _OPENMP  #endif  // _OPENMP
397  }  }
 #endif   // LAZY_NODE_STORAGE  
398    
399    
400  // Creates an identity node  // Creates an identity node
401  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
402      : parent(p->getFunctionSpace(),p->getShape())          : parent(p->getFunctionSpace(),p->getShape())
403  #ifdef LAZY_NODE_STORAGE          ,m_sampleids(0),
404      ,m_sampleids(0),          m_samples(1)
     m_samples(1)  
 #endif  
405  {  {
406     if (p->isLazy())     if (p->isLazy())
407     {     {
408      // I don't want identity of Lazy.          // I don't want identity of Lazy.
409      // Question: Why would that be so bad?          // Question: Why would that be so bad?
410      // Answer: We assume that the child of ID is something we can call getVector on          // Answer: We assume that the child of ID is something we can call getVector on
411      throw DataException("Programmer error - attempt to create identity from a DataLazy.");          throw DataException("Programmer error - attempt to create identity from a DataLazy.");
412     }     }
413     else     else
414     {     {
415      p->makeLazyShared();          DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
416      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);          makeIdentity(dr);
     makeIdentity(dr);  
417  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
418     }     }
419  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
420  }  }
421    
422  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
423      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),          : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
424      m_op(op),          m_op(op),
425      m_axis_offset(0),          m_axis_offset(0),
426      m_transpose(0),          m_transpose(0),
427      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
428  {  {
429     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
430     {     {
431      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");          throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
432     }     }
433    
434     DataLazy_ptr lleft;     DataLazy_ptr lleft;
435     if (!left->isLazy())     if (!left->isLazy())
436     {     {
437      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
438     }     }
439     else     else
440     {     {
441      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
442     }     }
443     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
444     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
445     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
446     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
447     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
448     LazyNodeSetup();     LazyNodeSetup();
 #endif  
449     SIZELIMIT     SIZELIMIT
450  }  }
451    
452    
453  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
454  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
455      : parent(resultFS(left,right,op), resultShape(left,right,op)),          : parent(resultFS(left,right,op), resultShape(left,right,op)),
456      m_op(op),          m_op(op),
457      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
458  {  {
459  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
460     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
461     {     {
462      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");          throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
463     }     }
464    
465     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
466     {     {
467      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
468      Data ltemp(left);          Data ltemp(left);
469      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
470      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
471     }     }
472     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
473     {     {
474      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
475      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
476  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
477     }     }
478     left->operandCheck(*right);     left->operandCheck(*right);
479    
480     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())                  // the children need to be DataLazy. Wrap them in IDENTITY if required
481     {     {
482      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
483  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
484     }     }
485     else     else
486     {     {
487      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
488  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
489     }     }
490     if (right->isLazy())     if (right->isLazy())
491     {     {
492      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
493  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
494     }     }
495     else     else
496     {     {
497      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
498  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
499     }     }
500     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
501     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
502     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
503     {     {
504      m_readytype='E';          m_readytype='E';
505     }     }
506     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
507     {     {
508      m_readytype='T';          m_readytype='T';
509     }     }
510     else     else
511     {     {
512      m_readytype='C';          m_readytype='C';
513     }     }
514     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());    
    m_buffsRequired=calcBuffs(m_left, m_right,m_op);  
515     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
516     m_height=max(m_left->m_height,m_right->m_height)+1;     m_height=max(m_left->m_height,m_right->m_height)+1;
 #ifdef LAZY_NODE_STORAGE  
517     LazyNodeSetup();     LazyNodeSetup();
 #endif  
518     SIZELIMIT     SIZELIMIT
519  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
520  }  }
521    
522  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
523      : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),          : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
524      m_op(op),          m_op(op),
525      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
526      m_transpose(transpose)          m_transpose(transpose)
527  {  {
528     if ((getOpgroup(op)!=G_TENSORPROD))     if ((getOpgroup(op)!=G_TENSORPROD))
529     {     {
530      throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");          throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
531     }     }
532     if ((transpose>2) || (transpose<0))     if ((transpose>2) || (transpose<0))
533     {     {
534      throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");          throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
535     }     }
536     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
537     {     {
538      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
539      Data ltemp(left);          Data ltemp(left);
540      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
541      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
542     }     }
543     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
544     {     {
545      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
546      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
547     }     }
548  //    left->operandCheck(*right);  //    left->operandCheck(*right);
549    
550     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())                  // the children need to be DataLazy. Wrap them in IDENTITY if required
551     {     {
552      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
553     }     }
554     else     else
555     {     {
556      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
557     }     }
558     if (right->isLazy())     if (right->isLazy())
559     {     {
560      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
561     }     }
562     else     else
563     {     {
564      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
565     }     }
566     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
567     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
568     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
569     {     {
570      m_readytype='E';          m_readytype='E';
571     }     }
572     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
573     {     {
574      m_readytype='T';          m_readytype='T';
575     }     }
576     else     else
577     {     {
578      m_readytype='C';          m_readytype='C';
579     }     }
580     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());    
    m_buffsRequired=calcBuffs(m_left, m_right,m_op);  
581     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
582     m_height=max(m_left->m_height,m_right->m_height)+1;     m_height=max(m_left->m_height,m_right->m_height)+1;
 #ifdef LAZY_NODE_STORAGE  
583     LazyNodeSetup();     LazyNodeSetup();
 #endif  
584     SIZELIMIT     SIZELIMIT
585  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
586  }  }
587    
588    
589  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
590      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),          : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
591      m_op(op),          m_op(op),
592      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
593      m_transpose(0),          m_transpose(0),
594      m_tol(0)          m_tol(0)
595  {  {
596     if ((getOpgroup(op)!=G_NP1OUT_P))     if ((getOpgroup(op)!=G_NP1OUT_P))
597     {     {
598      throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");          throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
599     }     }
600     DataLazy_ptr lleft;     DataLazy_ptr lleft;
601     if (!left->isLazy())     if (!left->isLazy())
602     {     {
603      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
604     }     }
605     else     else
606     {     {
607      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
608     }     }
609     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
610     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
611     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
612     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
613     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
614     LazyNodeSetup();     LazyNodeSetup();
 #endif  
615     SIZELIMIT     SIZELIMIT
616  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
617  }  }
618    
619  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
620      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
621      m_op(op),          m_op(op),
622      m_axis_offset(0),          m_axis_offset(0),
623      m_transpose(0),          m_transpose(0),
624      m_tol(tol)          m_tol(tol)
625  {  {
626     if ((getOpgroup(op)!=G_UNARY_P))     if ((getOpgroup(op)!=G_UNARY_P))
627     {     {
628      throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");          throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
629     }     }
630     DataLazy_ptr lleft;     DataLazy_ptr lleft;
631     if (!left->isLazy())     if (!left->isLazy())
632     {     {
633      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
634     }     }
635     else     else
636     {     {
637      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
638     }     }
639     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
640     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
641     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
642     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
643     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
644     LazyNodeSetup();     LazyNodeSetup();
 #endif  
645     SIZELIMIT     SIZELIMIT
646  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
647  }  }
648    
649    
650  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
651      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),          : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
652      m_op(op),          m_op(op),
653      m_axis_offset(axis0),          m_axis_offset(axis0),
654      m_transpose(axis1),          m_transpose(axis1),
655      m_tol(0)          m_tol(0)
656  {  {
657     if ((getOpgroup(op)!=G_NP1OUT_2P))     if ((getOpgroup(op)!=G_NP1OUT_2P))
658     {     {
659      throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");          throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
660     }     }
661     DataLazy_ptr lleft;     DataLazy_ptr lleft;
662     if (!left->isLazy())     if (!left->isLazy())
663     {     {
664      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
665     }     }
666     else     else
667     {     {
668      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
669     }     }
670     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
671     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
672     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
673     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
674     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
675     LazyNodeSetup();     LazyNodeSetup();
 #endif  
676     SIZELIMIT     SIZELIMIT
677  LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
678  }  }
679    
680  DataLazy::~DataLazy()  
681    namespace
682  {  {
 #ifdef LAZY_NODE_SETUP  
    delete[] m_sampleids;  
    delete[] m_samples;  
 #endif  
 }  
683    
684        inline int max3(int a, int b, int c)
685        {
686            int t=(a>b?a:b);
687            return (t>c?t:c);
688    
689  int      }
 DataLazy::getBuffsRequired() const  
 {  
     return m_buffsRequired;  
690  }  }
691    
692    DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
693  size_t          : parent(left->getFunctionSpace(), left->getShape()),
694  DataLazy::getMaxSampleSize() const          m_op(CONDEVAL),
695            m_axis_offset(0),
696            m_transpose(0),
697            m_tol(0)
698  {  {
699      return m_maxsamplesize;  
700       DataLazy_ptr lmask;
701       DataLazy_ptr lleft;
702       DataLazy_ptr lright;
703       if (!mask->isLazy())
704       {
705            lmask=DataLazy_ptr(new DataLazy(mask));
706       }
707       else
708       {
709            lmask=dynamic_pointer_cast<DataLazy>(mask);
710       }
711       if (!left->isLazy())
712       {
713            lleft=DataLazy_ptr(new DataLazy(left));
714       }
715       else
716       {
717            lleft=dynamic_pointer_cast<DataLazy>(left);
718       }
719       if (!right->isLazy())
720       {
721            lright=DataLazy_ptr(new DataLazy(right));
722       }
723       else
724       {
725            lright=dynamic_pointer_cast<DataLazy>(right);
726       }
727       m_readytype=lmask->m_readytype;
728       if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
729       {
730            throw DataException("Programmer Error - condEval arguments must have the same readytype");
731       }
732       m_left=lleft;
733       m_right=lright;
734       m_mask=lmask;
735       m_samplesize=getNumDPPSample()*getNoValues();
736       m_children=m_left->m_children+m_right->m_children+m_mask->m_children+1;
737       m_height=max3(m_left->m_height,m_right->m_height,m_mask->m_height)+1;
738       LazyNodeSetup();
739       SIZELIMIT
740    LAZYDEBUG(cout << "(8)Lazy created with " << m_samplesize << endl;)
741  }  }
742    
743    
744    
745  size_t  DataLazy::~DataLazy()
 DataLazy::getSampleBufferSize() const  
746  {  {
747      return m_maxsamplesize*(max(1,m_buffsRequired));     delete[] m_sampleids;
748  }  }
749    
750    
751  /*  /*
752    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
753    This does the work for the collapse method.    This does the work for the collapse method.
754    For reasons of efficiency do not call this method on DataExpanded nodes.    For reasons of efficiency do not call this method on DataExpanded nodes.
755  */  */
756  DataReady_ptr  DataReady_ptr
757  DataLazy::collapseToReady()  DataLazy::collapseToReady() const
758  {  {
759    if (m_readytype=='E')    if (m_readytype=='E')
760    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
761      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
762    }    }
763    if (m_op==IDENTITY)    if (m_op==IDENTITY)
# Line 823  DataLazy::collapseToReady() Line 775  DataLazy::collapseToReady()
775    switch(m_op)    switch(m_op)
776    {    {
777      case ADD:      case ADD:
778      result=left+right;          result=left+right;
779      break;          break;
780      case SUB:            case SUB:          
781      result=left-right;          result=left-right;
782      break;          break;
783      case MUL:            case MUL:          
784      result=left*right;          result=left*right;
785      break;          break;
786      case DIV:            case DIV:          
787      result=left/right;          result=left/right;
788      break;          break;
789      case SIN:      case SIN:
790      result=left.sin();            result=left.sin();      
791      break;          break;
792      case COS:      case COS:
793      result=left.cos();          result=left.cos();
794      break;          break;
795      case TAN:      case TAN:
796      result=left.tan();          result=left.tan();
797      break;          break;
798      case ASIN:      case ASIN:
799      result=left.asin();          result=left.asin();
800      break;          break;
801      case ACOS:      case ACOS:
802      result=left.acos();          result=left.acos();
803      break;          break;
804      case ATAN:      case ATAN:
805      result=left.atan();          result=left.atan();
806      break;          break;
807      case SINH:      case SINH:
808      result=left.sinh();          result=left.sinh();
809      break;          break;
810      case COSH:      case COSH:
811      result=left.cosh();          result=left.cosh();
812      break;          break;
813      case TANH:      case TANH:
814      result=left.tanh();          result=left.tanh();
815      break;          break;
816      case ERF:      case ERF:
817      result=left.erf();          result=left.erf();
818      break;          break;
819     case ASINH:     case ASINH:
820      result=left.asinh();          result=left.asinh();
821      break;          break;
822     case ACOSH:     case ACOSH:
823      result=left.acosh();          result=left.acosh();
824      break;          break;
825     case ATANH:     case ATANH:
826      result=left.atanh();          result=left.atanh();
827      break;          break;
828      case LOG10:      case LOG10:
829      result=left.log10();          result=left.log10();
830      break;          break;
831      case LOG:      case LOG:
832      result=left.log();          result=left.log();
833      break;          break;
834      case SIGN:      case SIGN:
835      result=left.sign();          result=left.sign();
836      break;          break;
837      case ABS:      case ABS:
838      result=left.abs();          result=left.abs();
839      break;          break;
840      case NEG:      case NEG:
841      result=left.neg();          result=left.neg();
842      break;          break;
843      case POS:      case POS:
844      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
845      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
846      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
847      break;          break;
848      case EXP:      case EXP:
849      result=left.exp();          result=left.exp();
850      break;          break;
851      case SQRT:      case SQRT:
852      result=left.sqrt();          result=left.sqrt();
853      break;          break;
854      case RECIP:      case RECIP:
855      result=left.oneOver();          result=left.oneOver();
856      break;          break;
857      case GZ:      case GZ:
858      result=left.wherePositive();          result=left.wherePositive();
859      break;          break;
860      case LZ:      case LZ:
861      result=left.whereNegative();          result=left.whereNegative();
862      break;          break;
863      case GEZ:      case GEZ:
864      result=left.whereNonNegative();          result=left.whereNonNegative();
865      break;          break;
866      case LEZ:      case LEZ:
867      result=left.whereNonPositive();          result=left.whereNonPositive();
868      break;          break;
869      case NEZ:      case NEZ:
870      result=left.whereNonZero(m_tol);          result=left.whereNonZero(m_tol);
871      break;          break;
872      case EZ:      case EZ:
873      result=left.whereZero(m_tol);          result=left.whereZero(m_tol);
874      break;          break;
875      case SYM:      case SYM:
876      result=left.symmetric();          result=left.symmetric();
877      break;          break;
878      case NSYM:      case NSYM:
879      result=left.nonsymmetric();          result=left.antisymmetric();
880      break;          break;
881      case PROD:      case PROD:
882      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);          result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
883      break;          break;
884      case TRANS:      case TRANS:
885      result=left.transpose(m_axis_offset);          result=left.transpose(m_axis_offset);
886      break;          break;
887      case TRACE:      case TRACE:
888      result=left.trace(m_axis_offset);          result=left.trace(m_axis_offset);
889      break;          break;
890      case SWAP:      case SWAP:
891      result=left.swapaxes(m_axis_offset, m_transpose);          result=left.swapaxes(m_axis_offset, m_transpose);
892      break;          break;
893      case MINVAL:      case MINVAL:
894      result=left.minval();          result=left.minval();
895      break;          break;
896      case MAXVAL:      case MAXVAL:
897      result=left.minval();          result=left.minval();
898            break;
899        case HER:
900        result=left.hermitian();
901      break;      break;
902      default:      default:
903      throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
904    }    }
905    return result.borrowReadyPtr();    return result.borrowReadyPtr();
906  }  }
# Line 957  DataLazy::collapseToReady() Line 912  DataLazy::collapseToReady()
912     the purpose of using DataLazy in the first place).     the purpose of using DataLazy in the first place).
913  */  */
914  void  void
915  DataLazy::collapse()  DataLazy::collapse() const
916  {  {
917    if (m_op==IDENTITY)    if (m_op==IDENTITY)
918    {    {
919      return;          return;
920    }    }
921    if (m_readytype=='E')    if (m_readytype=='E')
922    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
923      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
924    }    }
925    m_id=collapseToReady();    m_id=collapseToReady();
926    m_op=IDENTITY;    m_op=IDENTITY;
927  }  }
928    
 /*  
   \brief Compute the value of the expression (unary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
 {  
     // we assume that any collapsing has been done before we get here  
     // since we only have one argument we don't need to think about only  
     // processing single points.  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");  
   }  
   const ValueType* vleft=m_left->resolveVectorSample(v,offset,sampleNo,roffset);  
   const double* left=&((*vleft)[roffset]);  
   double* result=&(v[offset]);  
   roffset=offset;  
   switch (m_op)  
   {  
     case SIN:    
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);  
     break;  
     case COS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);  
     break;  
     case TAN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);  
     break;  
     case ASIN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);  
     break;  
     case ACOS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);  
     break;  
     case ATAN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);  
     break;  
     case SINH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);  
     break;  
     case COSH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);  
     break;  
     case TANH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);  
     break;  
     case ERF:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::erf);  
     break;  
 #endif  
    case ASINH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 #endif    
     break;  
    case ACOSH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 #endif    
     break;  
    case ATANH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 #endif    
     break;  
     case LOG10:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);  
     break;  
     case LOG:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);  
     break;  
     case SIGN:  
     tensor_unary_operation(m_samplesize, left, result, escript::fsign);  
     break;  
     case ABS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);  
     break;  
     case NEG:  
     tensor_unary_operation(m_samplesize, left, result, negate<double>());  
     break;  
     case POS:  
     // it doesn't mean anything for delayed.  
     // it will just trigger a deep copy of the lazy object  
     throw DataException("Programmer error - POS not supported for lazy data.");  
     break;  
     case EXP:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);  
     break;  
     case SQRT:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);  
     break;  
     case RECIP:  
     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));  
     break;  
     case GZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));  
     break;  
     case LZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));  
     break;  
     case GEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));  
     break;  
     case LEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));  
     break;  
 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently  
     case NEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));  
     break;  
     case EZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));  
     break;  
   
     default:  
     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");  
   }  
   return &v;  
 }  
   
   
 /*  
   \brief Compute the value of the expression (reduction operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveReduction(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
 {  
     // we assume that any collapsing has been done before we get here  
     // since we only have one argument we don't need to think about only  
     // processing single points.  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");  
   }  
   const ValueType* vleft=m_left->resolveVectorSample(v,offset,sampleNo,roffset);  
   double* result=&(v[offset]);  
   roffset=offset;  
   switch (m_op)  
   {  
     case MINVAL:  
     {  
       FMin op;  
       *result=DataMaths::reductionOp(*vleft, m_left->getShape(), roffset, op, numeric_limits<double>::max());  
     }  
     break;  
     case MAXVAL:  
     {  
       FMax op;  
       *result=DataMaths::reductionOp(*vleft, m_left->getShape(), roffset, op, numeric_limits<double>::max()*-1);  
     }  
     break;  
     default:  
     throw DataException("Programmer error - resolveReduction can not resolve operator "+opToString(m_op)+".");  
   }  
   return &v;  
 }  
   
   
   
 /*  
   \brief Compute the value of the expression (unary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
 {  
     // we assume that any collapsing has been done before we get here  
     // since we only have one argument we don't need to think about only  
     // processing single points.  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");  
   }  
     // since we can't write the result over the input, we need a result offset further along  
   size_t subroffset=roffset+m_samplesize;  
 LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)  
   const ValueType* vleft=m_left->resolveVectorSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);  
   roffset=offset;  
   size_t loop=0;  
   size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;  
   size_t step=getNoValues();  
   switch (m_op)  
   {  
     case SYM:  
     for (loop=0;loop<numsteps;++loop)  
     {  
         DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);  
         subroffset+=step;  
         offset+=step;  
     }  
     break;  
     case NSYM:  
     for (loop=0;loop<numsteps;++loop)  
     {  
         DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);  
         subroffset+=step;  
         offset+=step;  
     }  
     break;  
     default:  
     throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");  
   }  
   return &v;  
 }  
   
 /*  
   \brief Compute the value of the expression (unary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
 {  
     // we assume that any collapsing has been done before we get here  
     // since we only have one argument we don't need to think about only  
     // processing single points.  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");  
   }  
     // since we can't write the result over the input, we need a result offset further along  
   size_t subroffset;  
   const ValueType* vleft=m_left->resolveVectorSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);  
 LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)  
 LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)  
   roffset=offset;  
   size_t loop=0;  
   size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;  
   size_t outstep=getNoValues();  
   size_t instep=m_left->getNoValues();  
 LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)  
   switch (m_op)  
   {  
     case TRACE:  
     for (loop=0;loop<numsteps;++loop)  
     {  
 size_t zz=sampleNo*getNumDPPSample()+loop;  
 if (zz==5800)  
 {  
 LAZYDEBUG(cerr << "point=" <<  zz<< endl;)  
 LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)  
 LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)  
 LAZYDEBUG(cerr << subroffset << endl;)  
 LAZYDEBUG(cerr << "output=" << offset << endl;)  
 }  
             DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);  
 if (zz==5800)  
 {  
 LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)  
 }  
         subroffset+=instep;  
         offset+=outstep;  
     }  
     break;  
     case TRANS:  
     for (loop=0;loop<numsteps;++loop)  
     {  
             DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);  
         subroffset+=instep;  
         offset+=outstep;  
     }  
     break;  
     default:  
     throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");  
   }  
   return &v;  
 }  
   
   
 /*  
   \brief Compute the value of the expression (unary operation with int params) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT_2P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
 {  
     // we assume that any collapsing has been done before we get here  
     // since we only have one argument we don't need to think about only  
     // processing single points.  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer error - resolveNP1OUT_2P should only be called on expanded Data.");  
   }  
     // since we can't write the result over the input, we need a result offset further along  
   size_t subroffset;  
   const ValueType* vleft=m_left->resolveVectorSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);  
 LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)  
 LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)  
   roffset=offset;  
   size_t loop=0;  
   size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;  
   size_t outstep=getNoValues();  
   size_t instep=m_left->getNoValues();  
 LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)  
   switch (m_op)  
   {  
     case SWAP:  
     for (loop=0;loop<numsteps;++loop)  
     {  
             DataMaths::swapaxes(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset, m_transpose);  
         subroffset+=instep;  
         offset+=outstep;  
     }  
     break;  
     default:  
     throw DataException("Programmer error - resolveNP1OUT2P can not resolve operator "+opToString(m_op)+".");  
   }  
   return &v;  
 }  
   
   
   
 #define PROC_OP(TYPE,X)                               \  
     for (int j=0;j<onumsteps;++j)\  
     {\  
       for (int i=0;i<numsteps;++i,resultp+=resultStep) \  
       { \  
 LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\  
 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\  
          tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \  
 LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \  
          lroffset+=leftstep; \  
          rroffset+=rightstep; \  
       }\  
       lroffset+=oleftstep;\  
       rroffset+=orightstep;\  
     }  
   
 /*  
   \brief Compute the value of the expression (binary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  
 // have already been collapsed to IDENTITY. So we must have at least one expanded child.  
 // If both children are expanded, then we can process them in a single operation (we treat  
 // the whole sample as one big datapoint.  
 // If one of the children is not expanded, then we need to treat each point in the sample  
 // individually.  
 // There is an additional complication when scalar operations are considered.  
 // For example, 2+Vector.  
 // In this case each double within the point is treated individually  
 DataTypes::ValueType*  
 DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  
 {  
 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  
   
   size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors  
     // first work out which of the children are expanded  
   bool leftExp=(m_left->m_readytype=='E');  
   bool rightExp=(m_right->m_readytype=='E');  
   if (!leftExp && !rightExp)  
   {  
     throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");  
   }  
   bool leftScalar=(m_left->getRank()==0);  
   bool rightScalar=(m_right->getRank()==0);  
   if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))  
   {  
     throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");  
   }  
   size_t leftsize=m_left->getNoValues();  
   size_t rightsize=m_right->getNoValues();  
   size_t chunksize=1;           // how many doubles will be processed in one go  
   int leftstep=0;       // how far should the left offset advance after each step  
   int rightstep=0;  
   int numsteps=0;       // total number of steps for the inner loop  
   int oleftstep=0;  // the o variables refer to the outer loop  
   int orightstep=0; // The outer loop is only required in cases where there is an extended scalar  
   int onumsteps=1;  
     
   bool LES=(leftExp && leftScalar); // Left is an expanded scalar  
   bool RES=(rightExp && rightScalar);  
   bool LS=(!leftExp && leftScalar); // left is a single scalar  
   bool RS=(!rightExp && rightScalar);  
   bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar  
   bool RN=(!rightExp && !rightScalar);  
   bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar  
   bool REN=(rightExp && !rightScalar);  
   
   if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars  
   {  
     chunksize=m_left->getNumDPPSample()*leftsize;  
     leftstep=0;  
     rightstep=0;  
     numsteps=1;  
   }  
   else if (LES || RES)  
   {  
     chunksize=1;  
     if (LES)        // left is an expanded scalar  
     {  
         if (RS)  
         {  
            leftstep=1;  
            rightstep=0;  
            numsteps=m_left->getNumDPPSample();  
         }  
         else        // RN or REN  
         {  
            leftstep=0;  
            oleftstep=1;  
            rightstep=1;  
            orightstep=(RN ? -(int)rightsize : 0);  
            numsteps=rightsize;  
            onumsteps=m_left->getNumDPPSample();  
         }  
     }  
     else        // right is an expanded scalar  
     {  
         if (LS)  
         {  
            rightstep=1;  
            leftstep=0;  
            numsteps=m_right->getNumDPPSample();  
         }  
         else  
         {  
            rightstep=0;  
            orightstep=1;  
            leftstep=1;  
            oleftstep=(LN ? -(int)leftsize : 0);  
            numsteps=leftsize;  
            onumsteps=m_right->getNumDPPSample();  
         }  
     }  
   }  
   else  // this leaves (LEN, RS), (LEN, RN) and their transposes  
   {  
     if (LEN)    // and Right will be a single value  
     {  
         chunksize=rightsize;  
         leftstep=rightsize;  
         rightstep=0;  
         numsteps=m_left->getNumDPPSample();  
         if (RS)  
         {  
            numsteps*=leftsize;  
         }  
     }  
     else    // REN  
     {  
         chunksize=leftsize;  
         rightstep=leftsize;  
         leftstep=0;  
         numsteps=m_right->getNumDPPSample();  
         if (LS)  
         {  
            numsteps*=rightsize;  
         }  
     }  
   }  
   
   int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0  
     // Get the values of sub-expressions  
   const ValueType* left=m_left->resolveVectorSample(v,offset+getMaxSampleSize(),sampleNo,lroffset); // see note on  
     // calcBufss for why we can't put offset as the 2nd param above  
   const ValueType* right=m_right->resolveVectorSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note  
     // the right child starts further along.  
 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  
 LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)  
 LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)  
 LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)  
 LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)  
 LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)  
 LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)  
   
   
   double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved  
   switch(m_op)  
   {  
     case ADD:  
         PROC_OP(NO_ARG,plus<double>());  
     break;  
     case SUB:  
     PROC_OP(NO_ARG,minus<double>());  
     break;  
     case MUL:  
     PROC_OP(NO_ARG,multiplies<double>());  
     break;  
     case DIV:  
     PROC_OP(NO_ARG,divides<double>());  
     break;  
     case POW:  
        PROC_OP(double (double,double),::pow);  
     break;  
     default:  
     throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");  
   }  
   roffset=offset;  
   return &v;  
 }  
   
   
   
 /*  
   \brief Compute the value of the expression (tensor product) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  
 // have already been collapsed to IDENTITY. So we must have at least one expanded child.  
 // unlike the other resolve helpers, we must treat these datapoints separately.  
 DataTypes::ValueType*  
 DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  
 {  
 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)  
   
   size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors  
     // first work out which of the children are expanded  
   bool leftExp=(m_left->m_readytype=='E');  
   bool rightExp=(m_right->m_readytype=='E');  
   int steps=getNumDPPSample();  
 /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);  
   int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/  
   int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method  
   int rightStep=(rightExp?m_right->getNoValues() : 0);  
   
   int resultStep=getNoValues();  
     // Get the values of sub-expressions (leave a gap of one sample for the result).  
   int gap=offset+m_samplesize;    
   
 LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)  
   
   const ValueType* left=m_left->resolveVectorSample(v,gap,sampleNo,lroffset);  
   gap+=m_left->getMaxSampleSize();  
   
   
 LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)  
   
   
   const ValueType* right=m_right->resolveVectorSample(v,gap,sampleNo,rroffset);  
   
 LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;  
 cout << getNoValues() << endl;)  
 LAZYDEBUG(cerr << "Result of left=";)  
 LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;  
   
 for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)  
 {  
 cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";  
 if (i%4==0) cout << endl;  
 })  
 LAZYDEBUG(cerr << "\nResult of right=" << endl;)  
 LAZYDEBUG(  
 for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)  
 {  
 cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";  
 if (i%4==0) cout << endl;  
 }  
 cerr << endl;  
 )  
 LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)  
 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)  
 LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)  
 LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)  
 LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)  
 LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  
 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  
   
   double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved  
   switch(m_op)  
   {  
     case PROD:  
     for (int i=0;i<steps;++i,resultp+=resultStep)  
     {  
   
 LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)  
 LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)  
 LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)  
   
           const double *ptr_0 = &((*left)[lroffset]);  
           const double *ptr_1 = &((*right)[rroffset]);  
   
 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  
 LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  
   
           matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);  
   
 LAZYDEBUG(cout << "Results--\n";  
 {  
   DataVector dv(getNoValues());  
 for (int z=0;z<getNoValues();++z)  
 {  
   cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";  
   if (z%4==0) cout << endl;  
   dv[z]=resultp[z];  
 }  
 cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");  
 cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;  
 }  
 )  
       lroffset+=leftStep;  
       rroffset+=rightStep;  
     }  
     break;  
     default:  
     throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");  
   }  
   roffset=offset;  
   return &v;  
 }  
   
   
 #ifdef LAZY_NODE_STORAGE  
   
929  // The result will be stored in m_samples  // The result will be stored in m_samples
930  // The return value is a pointer to the DataVector, offset is the offset within the return value  // The return value is a pointer to the DataVector, offset is the offset within the return value
931  const DataTypes::ValueType*  const DataTypes::RealVectorType*
932  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
933  {  {
934  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
935      // collapse so we have a 'E' node or an IDENTITY for some other type          // collapse so we have a 'E' node or an IDENTITY for some other type
936    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
937    {    {
938      collapse();          collapse();
939    }    }
940    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
941    {    {
942      const ValueType& vec=m_id->getVectorRO();      const RealVectorType& vec=m_id->getVectorRO();
 //     if (m_readytype=='C')  
 //     {  
 //  roffset=0;      // all samples read from the same position  
 //  return &(m_samples);  
 //     }  
943      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
944    #ifdef LAZY_STACK_PROF
945    int x;
946    if (&x<stackend[omp_get_thread_num()])
947    {
948           stackend[omp_get_thread_num()]=&x;
949    }
950    #endif
951      return &(vec);      return &(vec);
952    }    }
953    if (m_readytype!='E')    if (m_readytype!='E')
# Line 1669  LAZYDEBUG(cout << "Resolve sample " << t Line 956  LAZYDEBUG(cout << "Resolve sample " << t
956    }    }
957    if (m_sampleids[tid]==sampleNo)    if (m_sampleids[tid]==sampleNo)
958    {    {
959      roffset=tid*m_samplesize;          roffset=tid*m_samplesize;
960      return &(m_samples);        // sample is already resolved          return &(m_samples);            // sample is already resolved
961    }    }
962    m_sampleids[tid]=sampleNo;    m_sampleids[tid]=sampleNo;
963    
964    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
965    {    {
966    case G_UNARY:    case G_UNARY:
# Line 1683  LAZYDEBUG(cout << "Resolve sample " << t Line 971  LAZYDEBUG(cout << "Resolve sample " << t
971    case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);    case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
972    case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);    case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
973    case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);    case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
974      case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
975    default:    default:
976      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
977    }    }
978  }  }
979    
980  const DataTypes::ValueType*  const DataTypes::RealVectorType*
981  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
982  {  {
983      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
984      // since we only have one argument we don't need to think about only          // since we only have one argument we don't need to think about only
985      // processing single points.          // processing single points.
986      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
987    if (m_readytype!='E')    if (m_readytype!='E')
988    {    {
989      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
# Line 1703  DataLazy::resolveNodeUnary(int tid, int Line 992  DataLazy::resolveNodeUnary(int tid, int
992    {    {
993      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
994    }    }
995    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
996    const double* left=&((*leftres)[roffset]);    const double* left=&((*leftres)[roffset]);
997    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
998    double* result=&(m_samples[roffset]);    double* result=&(m_samples[roffset]);
999    switch (m_op)    if (m_op==POS)
1000    {    {
1001      case SIN:        // this should be prevented earlier
1002      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);      // operation is meaningless for lazy
1003      break;          throw DataException("Programmer error - POS not supported for lazy data.");    
1004      case COS:    }
1005      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);    tensor_unary_array_operation(m_samplesize,
1006      break;                               left,
1007      case TAN:                               result,
1008      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);                               m_op,
1009      break;                               m_tol);  
     case ASIN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);  
     break;  
     case ACOS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);  
     break;  
     case ATAN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);  
     break;  
     case SINH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);  
     break;  
     case COSH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);  
     break;  
     case TANH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);  
     break;  
     case ERF:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::erf);  
     break;  
 #endif  
    case ASINH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 #endif    
     break;  
    case ACOSH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 #endif    
     break;  
    case ATANH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 #endif    
     break;  
     case LOG10:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);  
     break;  
     case LOG:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);  
     break;  
     case SIGN:  
     tensor_unary_operation(m_samplesize, left, result, escript::fsign);  
     break;  
     case ABS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);  
     break;  
     case NEG:  
     tensor_unary_operation(m_samplesize, left, result, negate<double>());  
     break;  
     case POS:  
     // it doesn't mean anything for delayed.  
     // it will just trigger a deep copy of the lazy object  
     throw DataException("Programmer error - POS not supported for lazy data.");  
     break;  
     case EXP:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);  
     break;  
     case SQRT:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);  
     break;  
     case RECIP:  
     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));  
     break;  
     case GZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));  
     break;  
     case LZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));  
     break;  
     case GEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));  
     break;  
     case LEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));  
     break;  
 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently  
     case NEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));  
     break;  
     case EZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));  
     break;  
   
     default:  
     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");  
   }  
1010    return &(m_samples);    return &(m_samples);
1011  }  }
1012    
1013    
1014  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1015  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1016  {  {
1017      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1018      // since we only have one argument we don't need to think about only          // since we only have one argument we don't need to think about only
1019      // processing single points.          // processing single points.
1020      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1021    if (m_readytype!='E')    if (m_readytype!='E')
1022    {    {
1023      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
# Line 1836  DataLazy::resolveNodeReduction(int tid, Line 1027  DataLazy::resolveNodeReduction(int tid,
1027      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1028    }    }
1029    size_t loffset=0;    size_t loffset=0;
1030    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1031    
1032    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1033      unsigned int ndpps=getNumDPPSample();
1034      unsigned int psize=DataTypes::noValues(m_left->getShape());
1035    double* result=&(m_samples[roffset]);    double* result=&(m_samples[roffset]);
1036    switch (m_op)    switch (m_op)
1037    {    {
1038      case MINVAL:      case MINVAL:
1039      {          {
1040        FMin op;            for (unsigned int z=0;z<ndpps;++z)
1041        *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());            {
1042      }              FMin op;
1043      break;              *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1044                loffset+=psize;
1045                result++;
1046              }
1047            }
1048            break;
1049      case MAXVAL:      case MAXVAL:
1050      {          {
1051        FMax op;            for (unsigned int z=0;z<ndpps;++z)
1052        *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);            {
1053      }            FMax op;
1054      break;            *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1055              loffset+=psize;
1056              result++;
1057              }
1058            }
1059            break;
1060      default:      default:
1061      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1062    }    }
1063    return &(m_samples);    return &(m_samples);
1064  }  }
1065    
1066  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1067  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1068  {  {
1069      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1070      // since we only have one argument we don't need to think about only          // since we only have one argument we don't need to think about only
1071      // processing single points.          // processing single points.
1072    if (m_readytype!='E')    if (m_readytype!='E')
1073    {    {
1074      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
# Line 1875  DataLazy::resolveNodeNP1OUT(int tid, int Line 1078  DataLazy::resolveNodeNP1OUT(int tid, int
1078      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1079    }    }
1080    size_t subroffset;    size_t subroffset;
1081    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1082    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1083    size_t loop=0;    size_t loop=0;
1084    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
# Line 1884  DataLazy::resolveNodeNP1OUT(int tid, int Line 1087  DataLazy::resolveNodeNP1OUT(int tid, int
1087    switch (m_op)    switch (m_op)
1088    {    {
1089      case SYM:      case SYM:
1090      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1091      {          {
1092          DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              escript::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1093          subroffset+=step;              subroffset+=step;
1094          offset+=step;              offset+=step;
1095      }          }
1096      break;          break;
1097      case NSYM:      case NSYM:
1098      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1099      {          {
1100          DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              escript::antisymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1101          subroffset+=step;              subroffset+=step;
1102          offset+=step;              offset+=step;
1103      }          }
1104      break;          break;
1105      default:      default:
1106      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1107    }    }
1108    return &m_samples;    return &m_samples;
1109  }  }
1110    
1111  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1112  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1113  {  {
1114      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1115      // since we only have one argument we don't need to think about only          // since we only have one argument we don't need to think about only
1116      // processing single points.          // processing single points.
1117    if (m_readytype!='E')    if (m_readytype!='E')
1118    {    {
1119      throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
# Line 1921  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1124  DataLazy::resolveNodeNP1OUT_P(int tid, i
1124    }    }
1125    size_t subroffset;    size_t subroffset;
1126    size_t offset;    size_t offset;
1127    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1128    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1129    offset=roffset;    offset=roffset;
1130    size_t loop=0;    size_t loop=0;
# Line 1931  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1134  DataLazy::resolveNodeNP1OUT_P(int tid, i
1134    switch (m_op)    switch (m_op)
1135    {    {
1136      case TRACE:      case TRACE:
1137      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1138      {          {
1139              DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);              escript::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1140          subroffset+=instep;              subroffset+=instep;
1141          offset+=outstep;              offset+=outstep;
1142      }          }
1143      break;          break;
1144      case TRANS:      case TRANS:
1145      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1146      {          {
1147              DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);              escript::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1148          subroffset+=instep;              subroffset+=instep;
1149          offset+=outstep;              offset+=outstep;
1150      }          }
1151      break;          break;
1152      default:      default:
1153      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1154    }    }
1155    return &m_samples;    return &m_samples;
1156  }  }
1157    
1158    
1159  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1160  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1161  {  {
1162    if (m_readytype!='E')    if (m_readytype!='E')
1163    {    {
# Line 1966  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1169  DataLazy::resolveNodeNP1OUT_2P(int tid,
1169    }    }
1170    size_t subroffset;    size_t subroffset;
1171    size_t offset;    size_t offset;
1172    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1173    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1174    offset=roffset;    offset=roffset;
1175    size_t loop=0;    size_t loop=0;
# Line 1976  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1179  DataLazy::resolveNodeNP1OUT_2P(int tid,
1179    switch (m_op)    switch (m_op)
1180    {    {
1181      case SWAP:      case SWAP:
1182      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1183      {          {
1184              DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);              escript::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1185          subroffset+=instep;              subroffset+=instep;
1186          offset+=outstep;              offset+=outstep;
1187      }          }
1188      break;          break;
1189      default:      default:
1190      throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1191    }    }
1192    return &m_samples;    return &m_samples;
1193  }  }
1194    
1195    const DataTypes::RealVectorType*
1196    DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1197    {
1198      if (m_readytype!='E')
1199      {
1200        throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1201      }
1202      if (m_op!=CONDEVAL)
1203      {
1204        throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1205      }
1206      size_t subroffset;
1207    
1208      const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1209      const RealVectorType* srcres=0;
1210      if ((*maskres)[subroffset]>0)
1211      {
1212            srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1213      }
1214      else
1215      {
1216            srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1217      }
1218    
1219      // Now we need to copy the result
1220    
1221      roffset=m_samplesize*tid;
1222      for (int i=0;i<m_samplesize;++i)
1223      {
1224            m_samples[roffset+i]=(*srcres)[subroffset+i];  
1225      }
1226    
1227      return &m_samples;
1228    }
1229    
1230  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1231  // have already been collapsed to IDENTITY. So we must have at least one expanded child.  // have already been collapsed to IDENTITY. So we must have at least one expanded child.
# Line 2000  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1236  DataLazy::resolveNodeNP1OUT_2P(int tid,
1236  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1237  // For example, 2+Vector.  // For example, 2+Vector.
1238  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1239  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1240  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1241  {  {
1242  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1243    
1244    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors    size_t lroffset=0, rroffset=0;        // offsets in the left and right result vectors
1245      // first work out which of the children are expanded          // first work out which of the children are expanded
1246    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1247    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1248    if (!leftExp && !rightExp)    if (!leftExp && !rightExp)
1249    {    {
1250      throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");          throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1251    }    }
1252    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1253    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1254    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1255    {    {
1256      throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1257    }    }
1258    size_t leftsize=m_left->getNoValues();    size_t leftsize=m_left->getNoValues();
1259    size_t rightsize=m_right->getNoValues();    size_t rightsize=m_right->getNoValues();
1260    size_t chunksize=1;           // how many doubles will be processed in one go    size_t chunksize=1;                   // how many doubles will be processed in one go
1261    int leftstep=0;       // how far should the left offset advance after each step    int leftstep=0;               // how far should the left offset advance after each step
1262    int rightstep=0;    int rightstep=0;
1263    int numsteps=0;       // total number of steps for the inner loop    int numsteps=0;               // total number of steps for the inner loop
1264    int oleftstep=0;  // the o variables refer to the outer loop    int oleftstep=0;      // the o variables refer to the outer loop
1265    int orightstep=0; // The outer loop is only required in cases where there is an extended scalar    int orightstep=0;     // The outer loop is only required in cases where there is an extended scalar
1266    int onumsteps=1;    int onumsteps=1;
1267        
1268    bool LES=(leftExp && leftScalar); // Left is an expanded scalar    bool LES=(leftExp && leftScalar);     // Left is an expanded scalar
1269    bool RES=(rightExp && rightScalar);    bool RES=(rightExp && rightScalar);
1270    bool LS=(!leftExp && leftScalar); // left is a single scalar    bool LS=(!leftExp && leftScalar);     // left is a single scalar
1271    bool RS=(!rightExp && rightScalar);    bool RS=(!rightExp && rightScalar);
1272    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1273    bool RN=(!rightExp && !rightScalar);    bool RN=(!rightExp && !rightScalar);
1274    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1275    bool REN=(rightExp && !rightScalar);    bool REN=(rightExp && !rightScalar);
1276    
1277    if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars    if ((LES && RES) || (LEN && REN))     // both are Expanded scalars or both are expanded non-scalars
1278    {    {
1279      chunksize=m_left->getNumDPPSample()*leftsize;          chunksize=m_left->getNumDPPSample()*leftsize;
1280      leftstep=0;          leftstep=0;
1281      rightstep=0;          rightstep=0;
1282      numsteps=1;          numsteps=1;
1283    }    }
1284    else if (LES || RES)    else if (LES || RES)
1285    {    {
1286      chunksize=1;          chunksize=1;
1287      if (LES)        // left is an expanded scalar          if (LES)                // left is an expanded scalar
1288      {          {
1289          if (RS)                  if (RS)
1290          {                  {
1291             leftstep=1;                     leftstep=1;
1292             rightstep=0;                     rightstep=0;
1293             numsteps=m_left->getNumDPPSample();                     numsteps=m_left->getNumDPPSample();
1294          }                  }
1295          else        // RN or REN                  else            // RN or REN
1296          {                  {
1297             leftstep=0;                     leftstep=0;
1298             oleftstep=1;                     oleftstep=1;
1299             rightstep=1;                     rightstep=1;
1300             orightstep=(RN ? -(int)rightsize : 0);                     orightstep=(RN ? -(int)rightsize : 0);
1301             numsteps=rightsize;                     numsteps=rightsize;
1302             onumsteps=m_left->getNumDPPSample();                     onumsteps=m_left->getNumDPPSample();
1303          }                  }
1304      }          }
1305      else        // right is an expanded scalar          else            // right is an expanded scalar
1306      {          {
1307          if (LS)                  if (LS)
1308          {                  {
1309             rightstep=1;                     rightstep=1;
1310             leftstep=0;                     leftstep=0;
1311             numsteps=m_right->getNumDPPSample();                     numsteps=m_right->getNumDPPSample();
1312          }                  }
1313          else                  else
1314          {                  {
1315             rightstep=0;                     rightstep=0;
1316             orightstep=1;                     orightstep=1;
1317             leftstep=1;                     leftstep=1;
1318             oleftstep=(LN ? -(int)leftsize : 0);                     oleftstep=(LN ? -(int)leftsize : 0);
1319             numsteps=leftsize;                     numsteps=leftsize;
1320             onumsteps=m_right->getNumDPPSample();                     onumsteps=m_right->getNumDPPSample();
1321          }                  }
1322      }          }
1323    }    }
1324    else  // this leaves (LEN, RS), (LEN, RN) and their transposes    else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1325    {    {
1326      if (LEN)    // and Right will be a single value          if (LEN)        // and Right will be a single value
1327      {          {
1328          chunksize=rightsize;                  chunksize=rightsize;
1329          leftstep=rightsize;                  leftstep=rightsize;
1330          rightstep=0;                  rightstep=0;
1331          numsteps=m_left->getNumDPPSample();                  numsteps=m_left->getNumDPPSample();
1332          if (RS)                  if (RS)
1333          {                  {
1334             numsteps*=leftsize;                     numsteps*=leftsize;
1335          }                  }
1336      }          }
1337      else    // REN          else    // REN
1338      {          {
1339          chunksize=leftsize;                  chunksize=leftsize;
1340          rightstep=leftsize;                  rightstep=leftsize;
1341          leftstep=0;                  leftstep=0;
1342          numsteps=m_right->getNumDPPSample();                  numsteps=m_right->getNumDPPSample();
1343          if (LS)                  if (LS)
1344          {                  {
1345             numsteps*=rightsize;                     numsteps*=rightsize;
1346          }                  }
1347      }          }
1348    }    }
1349    
1350    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1351      // Get the values of sub-expressions          // Get the values of sub-expressions
1352    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1353    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1354  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1355  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1356  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
# Line 2128  LAZYDEBUG(cout << "Right res["<< rroffse Line 1364  LAZYDEBUG(cout << "Right res["<< rroffse
1364    
1365    
1366    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1367    double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved    double* resultp=&(m_samples[roffset]);                // results are stored at the vector offset we received
1368    switch(m_op)    switch(m_op)
1369    {    {
1370      case ADD:      case ADD:
1371          PROC_OP(NO_ARG,plus<double>());          //PROC_OP(NO_ARG,plus<double>());
1372      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1373                 &(*left)[0],
1374                 &(*right)[0],
1375                 chunksize,
1376                 onumsteps,
1377                 numsteps,
1378                 resultStep,
1379                 leftstep,
1380                 rightstep,
1381                 oleftstep,
1382                 orightstep,
1383                 lroffset,
1384                 rroffset,
1385                 escript::ES_optype::ADD);  
1386            break;
1387      case SUB:      case SUB:
1388      PROC_OP(NO_ARG,minus<double>());        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1389      break;               &(*left)[0],
1390                 &(*right)[0],
1391                 chunksize,
1392                 onumsteps,
1393                 numsteps,
1394                 resultStep,
1395                 leftstep,
1396                 rightstep,
1397                 oleftstep,
1398                 orightstep,
1399                 lroffset,
1400                 rroffset,
1401                 escript::ES_optype::SUB);        
1402            //PROC_OP(NO_ARG,minus<double>());
1403            break;
1404      case MUL:      case MUL:
1405      PROC_OP(NO_ARG,multiplies<double>());          //PROC_OP(NO_ARG,multiplies<double>());
1406      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1407                 &(*left)[0],
1408                 &(*right)[0],
1409                 chunksize,
1410                 onumsteps,
1411                 numsteps,
1412                 resultStep,
1413                 leftstep,
1414                 rightstep,
1415                 oleftstep,
1416                 orightstep,
1417                 lroffset,
1418                 rroffset,
1419                 escript::ES_optype::MUL);        
1420            break;
1421      case DIV:      case DIV:
1422      PROC_OP(NO_ARG,divides<double>());          //PROC_OP(NO_ARG,divides<double>());
1423      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1424                 &(*left)[0],
1425                 &(*right)[0],
1426                 chunksize,
1427                 onumsteps,
1428                 numsteps,
1429                 resultStep,
1430                 leftstep,
1431                 rightstep,
1432                 oleftstep,
1433                 orightstep,
1434                 lroffset,
1435                 rroffset,
1436                 escript::ES_optype::DIV);        
1437            break;
1438      case POW:      case POW:
1439         PROC_OP(double (double,double),::pow);         //PROC_OP(double (double,double),::pow);
1440      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1441                 &(*left)[0],
1442                 &(*right)[0],
1443                 chunksize,
1444                 onumsteps,
1445                 numsteps,
1446                 resultStep,
1447                 leftstep,
1448                 rightstep,
1449                 oleftstep,
1450                 orightstep,
1451                 lroffset,
1452                 rroffset,
1453                 escript::ES_optype::POW);        
1454            break;
1455      default:      default:
1456      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1457    }    }
1458  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1459    return &m_samples;    return &m_samples;
# Line 2157  LAZYDEBUG(cout << "Result res[" << roffs Line 1463  LAZYDEBUG(cout << "Result res[" << roffs
1463  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1464  // have already been collapsed to IDENTITY. So we must have at least one expanded child.  // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1465  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1466  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1467  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1468  {  {
1469  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1470    
1471    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors    size_t lroffset=0, rroffset=0;        // offsets in the left and right result vectors
1472      // first work out which of the children are expanded          // first work out which of the children are expanded
1473    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1474    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1475    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1476    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method    int leftStep=(leftExp? m_left->getNoValues() : 0);            // do not have scalars as input to this method
1477    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1478    
1479    int resultStep=getNoValues();    int resultStep=getNoValues();
1480    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1481    size_t offset=roffset;    size_t offset=roffset;
1482    
1483    const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);    const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1484    
1485    const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1486    
1487  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1488  cout << getNoValues() << endl;)  cout << getNoValues() << endl;)
# Line 2190  LAZYDEBUG(cout << "m_samplesize=" << m_s Line 1496  LAZYDEBUG(cout << "m_samplesize=" << m_s
1496  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1497  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1498    
1499    double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(m_samples[offset]);         // results are stored at the vector offset we received
1500    switch(m_op)    switch(m_op)
1501    {    {
1502      case PROD:      case PROD:
1503      for (int i=0;i<steps;++i,resultp+=resultStep)          for (int i=0;i<steps;++i,resultp+=resultStep)
1504      {          {
1505            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1506            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1507    
1508  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1509  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1510    
1511            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1512    
1513        lroffset+=leftStep;            lroffset+=leftStep;
1514        rroffset+=rightStep;            rroffset+=rightStep;
1515      }          }
1516      break;          break;
1517      default:      default:
1518      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1519    }    }
1520    roffset=offset;    roffset=offset;
1521    return &m_samples;    return &m_samples;
1522  }  }
 #endif //LAZY_NODE_STORAGE  
   
 /*  
   \brief Compute the value of the expression for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
 */  
 // the vector and the offset are a place where the method could write its data if it wishes  
 // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
 // Hence the return value to indicate where the data is actually stored.  
 // Regardless, the storage should be assumed to be used, even if it isn't.  
   
 // the roffset is the offset within the returned vector where the data begins  
 const DataTypes::ValueType*  
 DataLazy::resolveVectorSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  
 {  
 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  
     // collapse so we have a 'E' node or an IDENTITY for some other type  
   if (m_readytype!='E' && m_op!=IDENTITY)  
   {  
     collapse();  
   }  
   if (m_op==IDENTITY)    
   {  
     const ValueType& vec=m_id->getVectorRO();  
     if (m_readytype=='C')  
     {  
     roffset=0;  
 LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)  
     return &(vec);  
     }  
     roffset=m_id->getPointOffset(sampleNo, 0);  
 LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)  
     return &(vec);  
   }  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer Error - Collapse did not produce an expanded node.");  
   }  
   switch (getOpgroup(m_op))  
   {  
   case G_UNARY:  
   case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);  
   case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);  
   case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);  
   case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);  
   case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);  
   case G_NP1OUT_2P: return resolveNP1OUT_2P(v, offset, sampleNo, roffset);  
   case G_REDUCTION: return resolveReduction(v, offset, sampleNo, roffset);  
   default:  
     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
   }  
1523    
 }  
1524    
1525  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1526  DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)  DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1527  {  {
1528  #ifdef _OPENMP  #ifdef _OPENMP
1529      int tid=omp_get_thread_num();          int tid=omp_get_thread_num();
1530  #else  #else
1531      int tid=0;          int tid=0;
1532  #endif  #endif
1533  #ifdef LAZY_NODE_STORAGE  
1534      return resolveNodeSample(tid, sampleNo, roffset);  #ifdef LAZY_STACK_PROF
1535            stackstart[tid]=&tid;
1536            stackend[tid]=&tid;
1537            const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1538            size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1539            #pragma omp critical
1540            if (d>maxstackuse)
1541            {
1542    cout << "Max resolve Stack use " << d << endl;
1543                    maxstackuse=d;
1544            }
1545            return r;
1546  #else  #else
1547      return resolveVectorSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);          return resolveNodeSample(tid, sampleNo, roffset);
1548  #endif  #endif
1549  }  }
1550    
# Line 2295  void Line 1554  void
1554  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
1555  {  {
1556     if (m_op==IDENTITY)     if (m_op==IDENTITY)
1557      return;          return;
 #ifndef LAZY_NODE_STORAGE  
    DataReady_ptr p=resolveVectorWorker();  
 #else  
1558     DataReady_ptr p=resolveNodeWorker();     DataReady_ptr p=resolveNodeWorker();
 #endif  
1559     makeIdentity(p);     makeIdentity(p);
1560  }  }
1561    
# Line 2316  void DataLazy::makeIdentity(const DataRe Line 1571  void DataLazy::makeIdentity(const DataRe
1571     else if(p->isExpanded()) {m_readytype='E';}     else if(p->isExpanded()) {m_readytype='E';}
1572     else if (p->isTagged()) {m_readytype='T';}     else if (p->isTagged()) {m_readytype='T';}
1573     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
    m_buffsRequired=1;  
1574     m_samplesize=p->getNumDPPSample()*p->getNoValues();     m_samplesize=p->getNumDPPSample()*p->getNoValues();
    m_maxsamplesize=m_samplesize;  
1575     m_left.reset();     m_left.reset();
1576     m_right.reset();     m_right.reset();
1577  }  }
# Line 2331  DataLazy::resolve() Line 1584  DataLazy::resolve()
1584      return m_id;      return m_id;
1585  }  }
1586    
 #ifdef LAZY_NODE_STORAGE  
1587    
1588  // This version of resolve uses storage in each node to hold results  /* This is really a static method but I think that caused problems in windows */
1589  DataReady_ptr  void
1590  DataLazy::resolveNodeWorker()  DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
1591  {  {
1592    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally    if (dats.empty())
1593    {    {
1594      collapse();          return;
   }  
   if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.  
   {  
     return m_id;  
1595    }    }
1596      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'    vector<DataLazy*> work;
1597    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    FunctionSpace fs=dats[0]->getFunctionSpace();
1598    ValueType& resvec=result->getVectorRW();    bool match=true;
1599    DataReady_ptr resptr=DataReady_ptr(result);    for (int i=dats.size()-1;i>=0;--i)
1600      {
1601    int sample;          if (dats[i]->m_readytype!='E')
1602    int totalsamples=getNumSamples();          {
1603    const ValueType* res=0;   // Storage for answer                  dats[i]->collapse();
1604  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)          }
1605    #pragma omp parallel for private(sample,res) schedule(static)          if (dats[i]->m_op!=IDENTITY)
1606    for (sample=0;sample<totalsamples;++sample)          {
1607    {                  work.push_back(dats[i]);
1608      size_t roffset=0;                  if (fs!=dats[i]->getFunctionSpace())
1609                    {
1610                            match=false;
1611                    }
1612            }
1613      }
1614      if (work.empty())
1615      {
1616            return;         // no work to do
1617      }
1618      if (match)    // all functionspaces match.  Yes I realise this is overly strict
1619      {             // it is possible that dats[0] is one of the objects which we discarded and
1620                    // all the other functionspaces match.
1621            vector<DataExpanded*> dep;
1622            vector<RealVectorType*> vecs;
1623            for (int i=0;i<work.size();++i)
1624            {
1625                    dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
1626                    vecs.push_back(&(dep[i]->getVectorRW()));
1627            }
1628            int totalsamples=work[0]->getNumSamples();
1629            const RealVectorType* res=0; // Storage for answer
1630            int sample;
1631            #pragma omp parallel private(sample, res)
1632            {
1633                size_t roffset=0;
1634                #pragma omp for schedule(static)
1635                for (sample=0;sample<totalsamples;++sample)
1636                {
1637                    roffset=0;
1638                    int j;
1639                    for (j=work.size()-1;j>=0;--j)
1640                    {
1641  #ifdef _OPENMP  #ifdef _OPENMP
1642      res=resolveNodeSample(omp_get_thread_num(),sample,roffset);                      res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1643  #else  #else
1644      res=resolveNodeSample(0,sample,roffset);                      res=work[j]->resolveNodeSample(0,sample,roffset);
1645  #endif  #endif
1646  LAZYDEBUG(cout << "Sample #" << sample << endl;)                      RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1647  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )                      memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
1648      DataVector::size_type outoffset=result->getPointOffset(sample,0);                  }
1649      memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));              }
1650            }
1651            // Now we need to load the new results as identity ops into the lazy nodes
1652            for (int i=work.size()-1;i>=0;--i)
1653            {
1654                work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1655            }
1656      }
1657      else  // functionspaces do not match
1658      {
1659            for (int i=0;i<work.size();++i)
1660            {
1661                    work[i]->resolveToIdentity();
1662            }
1663    }    }
   return resptr;  
1664  }  }
1665    
 #endif // LAZY_NODE_STORAGE  
1666    
1667  // To simplify the memory management, all threads operate on one large vector, rather than one each.  
1668  // Each sample is evaluated independently and copied into the result DataExpanded.  // This version of resolve uses storage in each node to hold results
1669  DataReady_ptr  DataReady_ptr
1670  DataLazy::resolveVectorWorker()  DataLazy::resolveNodeWorker()
1671  {  {
1672      if (m_readytype!='E')         // if the whole sub-expression is Constant or Tagged, then evaluate it normally
 LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)  
 LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)  
   if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally  
1673    {    {
1674      collapse();      collapse();
1675    }    }
1676    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.    if (m_op==IDENTITY)           // So a lazy expression of Constant or Tagged data will be returned here.
1677    {    {
1678      return m_id;      return m_id;
1679    }    }
1680      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'          // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1681    size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  RealVectorType(getNoValues()));
1682      // storage to evaluate its expression    RealVectorType& resvec=result->getVectorRW();
   int numthreads=1;  
 #ifdef _OPENMP  
   numthreads=omp_get_max_threads();  
 #endif  
   ValueType v(numthreads*threadbuffersize);  
 LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)  
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));  
   ValueType& resvec=result->getVectorRW();  
1683    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1684    
1685    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1686    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1687    const ValueType* res=0;   // Vector storing the answer    const RealVectorType* res=0;       // Storage for answer
   size_t resoffset=0;       // where in the vector to find the answer  
1688  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1689    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel private(sample,res)
   for (sample=0;sample<totalsamples;++sample)  
1690    {    {
1691  LAZYDEBUG(cout << "################################# " << sample << endl;)          size_t roffset=0;
1692    #ifdef LAZY_STACK_PROF
1693            stackstart[omp_get_thread_num()]=&roffset;
1694            stackend[omp_get_thread_num()]=&roffset;
1695    #endif
1696            #pragma omp for schedule(static)
1697            for (sample=0;sample<totalsamples;++sample)
1698            {
1699                    roffset=0;
1700  #ifdef _OPENMP  #ifdef _OPENMP
1701      res=resolveVectorSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);                  res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1702  #else  #else
1703      res=resolveVectorSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.                  res=resolveNodeSample(0,sample,roffset);
1704  #endif  #endif
1705  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1706  LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1707      outoffset=result->getPointOffset(sample,0);                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1708  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)                  memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1709      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector          }
1710      {    }
1711  LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)  #ifdef LAZY_STACK_PROF
1712      resvec[outoffset]=(*res)[resoffset];    for (int i=0;i<getNumberOfThreads();++i)
1713      }    {
1714  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)          size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1715  LAZYDEBUG(cerr << "*********************************" << endl;)  //      cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1716            if (r>maxstackuse)
1717            {
1718                    maxstackuse=r;
1719            }
1720    }    }
1721      cout << "Max resolve Stack use=" << maxstackuse << endl;
1722    #endif
1723    return resptr;    return resptr;
1724  }  }
1725    
# Line 2435  std::string Line 1727  std::string
1727  DataLazy::toString() const  DataLazy::toString() const
1728  {  {
1729    ostringstream oss;    ostringstream oss;
1730    oss << "Lazy Data:";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1731    intoString(oss);    switch (escriptParams.getLAZY_STR_FMT())
1732      {
1733      case 1:       // tree format
1734            oss << endl;
1735            intoTreeString(oss,"");
1736            break;
1737      case 2:       // just the depth
1738            break;
1739      default:
1740            intoString(oss);
1741            break;
1742      }
1743    return oss.str();    return oss.str();
1744  }  }
1745    
# Line 2448  DataLazy::intoString(ostringstream& oss) Line 1751  DataLazy::intoString(ostringstream& oss)
1751    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1752    {    {
1753    case G_IDENTITY:    case G_IDENTITY:
1754      if (m_id->isExpanded())          if (m_id->isExpanded())
1755      {          {
1756         oss << "E";             oss << "E";
1757      }          }
1758      else if (m_id->isTagged())          else if (m_id->isTagged())
1759      {          {
1760        oss << "T";            oss << "T";
1761      }          }
1762      else if (m_id->isConstant())          else if (m_id->isConstant())
1763      {          {
1764        oss << "C";            oss << "C";
1765      }          }
1766      else          else
1767      {          {
1768        oss << "?";            oss << "?";
1769      }          }
1770      oss << '@' << m_id.get();          oss << '@' << m_id.get();
1771      break;          break;
1772    case G_BINARY:    case G_BINARY:
1773      oss << '(';          oss << '(';
1774      m_left->intoString(oss);          m_left->intoString(oss);
1775      oss << ' ' << opToString(m_op) << ' ';          oss << ' ' << opToString(m_op) << ' ';
1776      m_right->intoString(oss);          m_right->intoString(oss);
1777      oss << ')';          oss << ')';
1778      break;          break;
1779    case G_UNARY:    case G_UNARY:
1780    case G_UNARY_P:    case G_UNARY_P:
1781    case G_NP1OUT:    case G_NP1OUT:
1782    case G_NP1OUT_P:    case G_NP1OUT_P:
1783    case G_REDUCTION:    case G_REDUCTION:
1784      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1785      m_left->intoString(oss);          m_left->intoString(oss);
1786      oss << ')';          oss << ')';
1787      break;          break;
1788    case G_TENSORPROD:    case G_TENSORPROD:
1789      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1790      m_left->intoString(oss);          m_left->intoString(oss);
1791      oss << ", ";          oss << ", ";
1792      m_right->intoString(oss);          m_right->intoString(oss);
1793      oss << ')';          oss << ')';
1794      break;          break;
1795    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1796      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1797      m_left->intoString(oss);          m_left->intoString(oss);
1798      oss << ", " << m_axis_offset << ", " << m_transpose;          oss << ", " << m_axis_offset << ", " << m_transpose;
1799      oss << ')';          oss << ')';
1800      break;          break;
1801      case G_CONDEVAL:
1802            oss << opToString(m_op)<< '(' ;
1803            m_mask->intoString(oss);
1804            oss << " ? ";
1805            m_left->intoString(oss);
1806            oss << " : ";
1807            m_right->intoString(oss);
1808            oss << ')';
1809            break;
1810    default:    default:
1811      oss << "UNKNOWN";          oss << "UNKNOWN";
1812    }    }
1813  }  }
1814    
1815    
1816    void
1817    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1818    {
1819      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1820      switch (getOpgroup(m_op))
1821      {
1822      case G_IDENTITY:
1823            if (m_id->isExpanded())
1824            {
1825               oss << "E";
1826            }
1827            else if (m_id->isTagged())
1828            {
1829              oss << "T";
1830            }
1831            else if (m_id->isConstant())
1832            {
1833              oss << "C";
1834            }
1835            else
1836            {
1837              oss << "?";
1838            }
1839            oss << '@' << m_id.get() << endl;
1840            break;
1841      case G_BINARY:
1842            oss << opToString(m_op) << endl;
1843            indent+='.';
1844            m_left->intoTreeString(oss, indent);
1845            m_right->intoTreeString(oss, indent);
1846            break;
1847      case G_UNARY:
1848      case G_UNARY_P:
1849      case G_NP1OUT:
1850      case G_NP1OUT_P:
1851      case G_REDUCTION:
1852            oss << opToString(m_op) << endl;
1853            indent+='.';
1854            m_left->intoTreeString(oss, indent);
1855            break;
1856      case G_TENSORPROD:
1857            oss << opToString(m_op) << endl;
1858            indent+='.';
1859            m_left->intoTreeString(oss, indent);
1860            m_right->intoTreeString(oss, indent);
1861            break;
1862      case G_NP1OUT_2P:
1863            oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1864            indent+='.';
1865            m_left->intoTreeString(oss, indent);
1866            break;
1867      default:
1868            oss << "UNKNOWN";
1869      }
1870    }
1871    
1872    
1873  DataAbstract*  DataAbstract*
1874  DataLazy::deepCopy()  DataLazy::deepCopy() const
1875  {  {
1876    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1877    {    {
1878    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1879    case G_UNARY:    case G_UNARY:
1880    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1881    case G_UNARY_P:   return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);    case G_UNARY_P:       return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1882    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);    case G_BINARY:        return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1883    case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);    case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1884    case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);    case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1885    case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);    case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
1886    case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);    case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1887    default:    default:
1888      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");          throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1889    }    }
1890  }  }
1891    
# Line 2527  DataLazy::deepCopy() Line 1897  DataLazy::deepCopy()
1897  // or it could be some function of the lengths of the DataReady instances which  // or it could be some function of the lengths of the DataReady instances which
1898  // form part of the expression.  // form part of the expression.
1899  // Rather than have people making assumptions, I have disabled the method.  // Rather than have people making assumptions, I have disabled the method.
1900  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1901  DataLazy::getLength() const  DataLazy::getLength() const
1902  {  {
1903    throw DataException("getLength() does not make sense for lazy data.");    throw DataException("getLength() does not make sense for lazy data.");
# Line 2542  DataLazy::getSlice(const DataTypes::Regi Line 1912  DataLazy::getSlice(const DataTypes::Regi
1912    
1913    
1914  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
1915  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1916  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1917                   int dataPointNo)                   int dataPointNo)
1918  {  {
1919    if (m_op==IDENTITY)    if (m_op==IDENTITY)
1920    {    {
1921      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1922    }    }
1923    if (m_readytype!='E')    if (m_readytype!='E')
1924    {    {
1925      collapse();          collapse();
1926      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1927    }    }
1928    // at this point we do not have an identity node and the expression will be Expanded    // at this point we do not have an identity node and the expression will be Expanded
1929    // so we only need to know which child to ask    // so we only need to know which child to ask
1930    if (m_left->m_readytype=='E')    if (m_left->m_readytype=='E')
1931    {    {
1932      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
1933    }    }
1934    else    else
1935    {    {
1936      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
1937    }    }
1938  }  }
1939    
1940  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
1941  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1942  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1943                   int dataPointNo) const                   int dataPointNo) const
1944  {  {
1945    if (m_op==IDENTITY)    if (m_op==IDENTITY)
1946    {    {
1947      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1948    }    }
1949    if (m_readytype=='E')    if (m_readytype=='E')
1950    {    {
# Line 2582  DataLazy::getPointOffset(int sampleNo, Line 1952  DataLazy::getPointOffset(int sampleNo,
1952      // so we only need to know which child to ask      // so we only need to know which child to ask
1953      if (m_left->m_readytype=='E')      if (m_left->m_readytype=='E')
1954      {      {
1955      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
1956      }      }
1957      else      else
1958      {      {
1959      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
1960      }      }
1961    }    }
1962    if (m_readytype=='C')    if (m_readytype=='C')
1963    {    {
1964      return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter          return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1965    }    }
1966    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1967  }  }
# Line 2601  DataLazy::getPointOffset(int sampleNo, Line 1971  DataLazy::getPointOffset(int sampleNo,
1971  void  void
1972  DataLazy::setToZero()  DataLazy::setToZero()
1973  {  {
1974  //   DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::RealVectorType v(getNoValues(),0);
1975  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1976  //   m_op=IDENTITY;  //   m_op=IDENTITY;
1977  //   m_right.reset();    //   m_right.reset();  
# Line 2609  DataLazy::setToZero() Line 1979  DataLazy::setToZero()
1979  //   m_readytype='C';  //   m_readytype='C';
1980  //   m_buffsRequired=1;  //   m_buffsRequired=1;
1981    
1982    privdebug=privdebug;  // to stop the compiler complaining about unused privdebug    (void)privdebug;  // to stop the compiler complaining about unused privdebug
1983    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1984  }  }
1985    
1986  bool  bool
1987  DataLazy::actsExpanded() const  DataLazy::actsExpanded() const
1988  {  {
1989      return (m_readytype=='E');          return (m_readytype=='E');
1990  }  }
1991    
1992  }   // end namespace  } // end namespace
1993    

Legend:
Removed from v.2721  
changed lines
  Added in v.6168

  ViewVC Help
Powered by ViewVC 1.1.26