/[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 2514 by jfenwick, Fri Jul 3 00:57:45 2009 UTC branches/clazy/escriptcore/src/DataLazy.cpp revision 6512 by jfenwick, Fri Mar 3 02:01:57 2017 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2008 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 42  bool privdebug=false; Line 38  bool privdebug=false;
38  #define DISABLEDEBUG privdebug=false;  #define DISABLEDEBUG privdebug=false;
39  }  }
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.getInt("TOO_MANY_LEVELS")) || (m_children>escript::escriptParams.getInt("TOO_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.getInt("TOO_MANY_LEVELS")) || (m_children>escript::escriptParams.getInt("TOO_MANY_NODES"))) {cerr << "SIZE LIMIT EXCEEDED " << m_height << endl;resolveToIdentity();}
44    
45    
46    #define SIZELIMIT \
47        if (m_height > escript::escriptParams.getTooManyLevels()) {\
48            if (escript::escriptParams.getLazyVerbose()) {\
49                cerr << "SIZE LIMIT EXCEEDED height=" << m_height << endl;\
50            }\
51            resolveToIdentity();\
52        }
53    
54  /*  /*
55  How does DataLazy work?  How does DataLazy work?
56  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 58  A special operation, IDENTITY, stores an Line 62  A special operation, IDENTITY, stores an
62  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.
63    
64  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, ...
65  Note that IDENITY is not considered a unary operation.  Note that IDENTITY is not considered a unary operation.
66    
67  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).
68  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 70  I will refer to individual DataLazy obje
70    
71  Each node also stores:  Each node also stores:
72  - 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
73      evaluated.          evaluated.
74  - 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
75      evaluate the expression.          evaluate the expression.
76  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
77    
78  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 113  namespace escript
113  namespace  namespace
114  {  {
115    
 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  
 };  
116    
117    // enabling this will print out when ever the maximum stacksize used by resolve increases
118    // it assumes _OPENMP is also in use
119    //#define LAZY_STACK_PROF
120    
121    
122    
123  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",  #ifndef _OPENMP
124              "sin","cos","tan",    #ifdef LAZY_STACK_PROF
125              "asin","acos","atan","sinh","cosh","tanh","erf",    #undef LAZY_STACK_PROF
126              "asinh","acosh","atanh",    #endif
127              "log10","log","sign","abs","neg","pos","exp","sqrt",  #endif
128              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",  
129              "symmetric","nonsymmetric",  
130              "prod",  #ifdef LAZY_STACK_PROF
131              "transpose", "trace",  std::vector<void*> stackstart(getNumberOfThreads());
132              "swapaxes"};  std::vector<void*> stackend(getNumberOfThreads());
133  int ES_opcount=41;  size_t maxstackuse=0;
134  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,  #endif
135              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};  
 inline  
 ES_opgroup  
 getOpgroup(ES_optype op)  
 {  
   return opgroups[op];  
 }  
136    
137  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
138  FunctionSpace  FunctionSpace
139  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
140  {  {
141      // perhaps this should call interpolate and throw or something?          // perhaps this should call interpolate and throw or something?
142      // maybe we need an interpolate node -          // maybe we need an interpolate node -
143      // that way, if interpolate is required in any other op we can just throw a          // that way, if interpolate is required in any other op we can just throw a
144      // programming error exception.          // programming error exception.
145    
146    FunctionSpace l=left->getFunctionSpace();    FunctionSpace l=left->getFunctionSpace();
147    FunctionSpace r=right->getFunctionSpace();    FunctionSpace r=right->getFunctionSpace();
148    if (l!=r)    if (l!=r)
149    {    {
150      if (r.probeInterpolation(l))      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
151        if (res==1)
152      {      {
153      return l;          return l;
154      }      }
155      if (l.probeInterpolation(r))      if (res==-1)
156      {      {
157      return r;          return r;
158      }      }
159      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
160    }    }
# Line 184  resultFS(DataAbstract_ptr left, DataAbst Line 166  resultFS(DataAbstract_ptr left, DataAbst
166  DataTypes::ShapeType  DataTypes::ShapeType
167  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
168  {  {
169      if (left->getShape()!=right->getShape())          if (left->getShape()!=right->getShape())
170      {          {
171        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))            if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
172        {            {
173          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.");
174        }            }
175        if (left->getRank()==0)   // we need to allow scalar * anything  
176        {            if (left->getRank()==0)       // we need to allow scalar * anything
177          return right->getShape();            {
178        }                  return right->getShape();
179        if (right->getRank()==0)            }
180        {            if (right->getRank()==0)
181          return left->getShape();            {
182        }                  return left->getShape();
183        throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");            }
184      }            throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
185      return left->getShape();          }
186            return left->getShape();
187  }  }
188    
189  // return the shape for "op left"  // return the shape for "op left"
# Line 208  resultShape(DataAbstract_ptr left, DataA Line 191  resultShape(DataAbstract_ptr left, DataA
191  DataTypes::ShapeType  DataTypes::ShapeType
192  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
193  {  {
194      switch(op)          switch(op)
195      {          {
196          case TRANS:          case TRANS:
197         {            // for the scoping of variables             {                    // for the scoping of variables
198          const DataTypes::ShapeType& s=left->getShape();                  const DataTypes::ShapeType& s=left->getShape();
199          DataTypes::ShapeType sh;                  DataTypes::ShapeType sh;
200          int rank=left->getRank();                  int rank=left->getRank();
201          if (axis_offset<0 || axis_offset>rank)                  if (axis_offset<0 || axis_offset>rank)
202          {                  {
203                 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);              stringstream e;
204              }              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
205              for (int i=0; i<rank; i++)              throw DataException(e.str());
206          {          }
207             int index = (axis_offset+i)%rank;          for (int i=0; i<rank; i++)
208                 sh.push_back(s[index]); // Append to new shape                  {
209              }                     int index = (axis_offset+i)%rank;
210          return sh;             sh.push_back(s[index]); // Append to new shape
211         }          }
212      break;                  return sh;
213      case TRACE:             }
214         {          break;
215          int rank=left->getRank();          case TRACE:
216          if (rank<2)             {
217          {                  int rank=left->getRank();
218             throw DataException("Trace can only be computed for objects with rank 2 or greater.");                  if (rank<2)
219          }                  {
220          if ((axis_offset>rank-2) || (axis_offset<0))                     throw DataException("Trace can only be computed for objects with rank 2 or greater.");
221          {                  }
222             throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");                  if ((axis_offset>rank-2) || (axis_offset<0))
223          }                  {
224          if (rank==2)                     throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
225          {                  }
226             return DataTypes::scalarShape;                  if (rank==2)
227          }                  {
228          else if (rank==3)                     return DataTypes::scalarShape;
229          {                  }
230             DataTypes::ShapeType sh;                  else if (rank==3)
231                 if (axis_offset==0)                  {
232             {                     DataTypes::ShapeType sh;
233                  sh.push_back(left->getShape()[2]);                     if (axis_offset==0)
234                 }                     {
235                 else     // offset==1                          sh.push_back(left->getShape()[2]);
236             {                     }
237              sh.push_back(left->getShape()[0]);                     else         // offset==1
238                 }                     {
239             return sh;                          sh.push_back(left->getShape()[0]);
240          }                     }
241          else if (rank==4)                     return sh;
242          {                  }
243             DataTypes::ShapeType sh;                  else if (rank==4)
244             const DataTypes::ShapeType& s=left->getShape();                  {
245                 if (axis_offset==0)                     DataTypes::ShapeType sh;
246             {                     const DataTypes::ShapeType& s=left->getShape();
247                  sh.push_back(s[2]);                     if (axis_offset==0)
248                  sh.push_back(s[3]);                     {
249                 }                          sh.push_back(s[2]);
250                 else if (axis_offset==1)                          sh.push_back(s[3]);
251             {                     }
252                  sh.push_back(s[0]);                     else if (axis_offset==1)
253                  sh.push_back(s[3]);                     {
254                 }                          sh.push_back(s[0]);
255             else     // offset==2                          sh.push_back(s[3]);
256             {                     }
257              sh.push_back(s[0]);                     else         // offset==2
258              sh.push_back(s[1]);                     {
259             }                          sh.push_back(s[0]);
260             return sh;                          sh.push_back(s[1]);
261          }                     }
262          else        // unknown rank                     return sh;
263          {                  }
264             throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");                  else            // unknown rank
265          }                  {
266         }                     throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
267      break;                  }
268          default:             }
269      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");          break;
270      }          default:
271            throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
272            }
273  }  }
274    
275  DataTypes::ShapeType  DataTypes::ShapeType
# Line 302  SwapShape(DataAbstract_ptr left, const i Line 287  SwapShape(DataAbstract_ptr left, const i
287          throw DataException("Error - Data::swapaxes argument must have at least rank 2.");          throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
288       }       }
289       if (axis0<0 || axis0>rank-1) {       if (axis0<0 || axis0>rank-1) {
290          throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);          stringstream e;
291            e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
292            throw DataException(e.str());
293       }       }
294       if (axis1<0 || axis1>rank-1) {       if (axis1<0 || axis1>rank-1) {
295           throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);          stringstream e;
296            e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
297            throw DataException(e.str());
298       }       }
299       if (axis0 == axis1) {       if (axis0 == axis1) {
300           throw DataException("Error - Data::swapaxes: axis indices must be different.");           throw DataException("Error - Data::swapaxes: axis indices must be different.");
# Line 336  SwapShape(DataAbstract_ptr left, const i Line 325  SwapShape(DataAbstract_ptr left, const i
325  DataTypes::ShapeType  DataTypes::ShapeType
326  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)
327  {  {
328                
329    // Get rank and shape of inputs    // Get rank and shape of inputs
330    int rank0 = left->getRank();    int rank0 = left->getRank();
331    int rank1 = right->getRank();    int rank1 = right->getRank();
# Line 345  GTPShape(DataAbstract_ptr left, DataAbst Line 334  GTPShape(DataAbstract_ptr left, DataAbst
334    
335    // Prepare for the loops of the product and verify compatibility of shapes    // Prepare for the loops of the product and verify compatibility of shapes
336    int start0=0, start1=0;    int start0=0, start1=0;
337    if (transpose == 0)       {}    if (transpose == 0)           {}
338    else if (transpose == 1)  { start0 = axis_offset; }    else if (transpose == 1)      { start0 = axis_offset; }
339    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)      { start1 = rank1-axis_offset; }
340    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"); }
341    
342    if (rank0<axis_offset)    if (rank0<axis_offset)
343    {    {
344      throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");          throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
345    }    }
346    
347    // Adjust the shapes for transpose    // Adjust the shapes for transpose
348    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather    DataTypes::ShapeType tmpShape0(rank0);        // pre-sizing the vectors rather
349    DataTypes::ShapeType tmpShape1(rank1);    // than using push_back    DataTypes::ShapeType tmpShape1(rank1);        // than using push_back
350    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]; }
351    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]; }
352    
353    // Prepare for the loops of the product    // Prepare for the loops of the product
354    SL=1, SM=1, SR=1;    SL=1, SM=1, SR=1;
355    for (int i=0; i<rank0-axis_offset; i++)   {    for (int i=0; i<rank0-axis_offset; i++)       {
356      SL *= tmpShape0[i];      SL *= tmpShape0[i];
357    }    }
358    for (int i=rank0-axis_offset; i<rank0; i++)   {    for (int i=rank0-axis_offset; i<rank0; i++)   {
359      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
360        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
361      }      }
362      SM *= tmpShape0[i];      SM *= tmpShape0[i];
363    }    }
364    for (int i=axis_offset; i<rank1; i++)     {    for (int i=axis_offset; i<rank1; i++)         {
365      SR *= tmpShape1[i];      SR *= tmpShape1[i];
366    }    }
367    
368    // 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)
369    DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      
370    {         // block to limit the scope of out_index    {                     // block to limit the scope of out_index
371       int out_index=0;       int out_index=0;
372       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
373       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 394  GTPShape(DataAbstract_ptr left, DataAbst Line 383  GTPShape(DataAbstract_ptr left, DataAbst
383    return shape2;    return shape2;
384  }  }
385    
386  // 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_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];  
 }  
387    
 #ifdef LAZY_NODE_STORAGE  
388  void DataLazy::LazyNodeSetup()  void DataLazy::LazyNodeSetup()
389  {  {
390  #ifdef _OPENMP  #ifdef _OPENMP
391      int numthreads=omp_get_max_threads();      int numthreads=omp_get_max_threads();
392      m_samples.resize(numthreads*m_samplesize);      m_samples_r.resize(numthreads*m_samplesize);
393      m_sampleids=new int[numthreads];      m_sampleids=new int[numthreads];
394      for (int i=0;i<numthreads;++i)      for (int i=0;i<numthreads;++i)
395      {      {
396          m_sampleids[i]=-1;            m_sampleids[i]=-1;  
397      }      }
398  #else  #else
399      m_samples.resize(m_samplesize);      if (m_iscompl)
400        {
401            m_samples_c.resize(m_samplesize);
402        }
403        else
404        {
405            m_samples_r.resize(m_samplesize);
406        }
407      m_sampleids=new int[1];      m_sampleids=new int[1];
408      m_sampleids[0]=-1;      m_sampleids[0]=-1;
409  #endif  // _OPENMP  #endif  // _OPENMP
410  }  }
 #endif   // LAZY_NODE_STORAGE  
411    
412    
413  // Creates an identity node  // Creates an identity node
414  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
415      : parent(p->getFunctionSpace(),p->getShape())          : parent(p->getFunctionSpace(),p->getShape())
416  #ifdef LAZY_NODE_STORAGE          ,m_sampleids(0),
417      ,m_sampleids(0),          m_samples_r(1)
     m_samples(1)  
 #endif  
418  {  {
419     if (p->isLazy())     if (p->isLazy())
420     {     {
421      // I don't want identity of Lazy.          // I don't want identity of Lazy.
422      // Question: Why would that be so bad?          // Question: Why would that be so bad?
423      // 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
424      throw DataException("Programmer error - attempt to create identity from a DataLazy.");          throw DataException("Programmer error - attempt to create identity from a DataLazy.");
425     }     }
426     else     else
427     {     {
428      p->makeLazyShared();          DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
429      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);          makeIdentity(dr);
     makeIdentity(dr);  
430  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
431     }     }
432  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
433  }  }
434    
435  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
436      : parent(left->getFunctionSpace(),left->getShape()),          : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
437      m_op(op),          m_op(op),
438      m_axis_offset(0),          m_axis_offset(0),
439      m_transpose(0),          m_transpose(0),
440      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
441  {  {
442     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
443     {     {
444      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.");
445     }     }
446    
447     DataLazy_ptr lleft;     DataLazy_ptr lleft;
448     if (!left->isLazy())     if (!left->isLazy())
449     {     {
450      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
451     }     }
452     else     else
453     {     {
454      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
455     }     }
456     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
457     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
458     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
459     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
460     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
461     LazyNodeSetup();     LazyNodeSetup();
 #endif  
462     SIZELIMIT     SIZELIMIT
463  }  }
464    
465    
466  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
467  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
468      : parent(resultFS(left,right,op), resultShape(left,right,op)),          : parent(resultFS(left,right,op), resultShape(left,right,op)),
469      m_op(op),          m_op(op),
470      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
471  {  {
472  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
473     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
474     {     {
475      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.");
476     }     }
477    
478     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
479     {     {
480      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
481      Data ltemp(left);          Data ltemp(left);
482      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
483      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
484     }     }
485     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
486     {     {
487      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
488      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
489  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
490     }     }
491     left->operandCheck(*right);     left->operandCheck(*right);
492    
493     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
494     {     {
495      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
496  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
497     }     }
498     else     else
499     {     {
500      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
501  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
502     }     }
503     if (right->isLazy())     if (right->isLazy())
504     {     {
505      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
506  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
507     }     }
508     else     else
509     {     {
510      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
511  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
512     }     }
513     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
514     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
515     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
516     {     {
517      m_readytype='E';          m_readytype='E';
518     }     }
519     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
520     {     {
521      m_readytype='T';          m_readytype='T';
522     }     }
523     else     else
524     {     {
525      m_readytype='C';          m_readytype='C';
526     }     }
527     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);  
528     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
529     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  
530     LazyNodeSetup();     LazyNodeSetup();
 #endif  
531     SIZELIMIT     SIZELIMIT
532  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
533  }  }
534    
535  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)
536      : 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)),
537      m_op(op),          m_op(op),
538      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
539      m_transpose(transpose)          m_transpose(transpose)
540  {  {
541     if ((getOpgroup(op)!=G_TENSORPROD))     if ((getOpgroup(op)!=G_TENSORPROD))
542     {     {
543      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.");
544     }     }
545     if ((transpose>2) || (transpose<0))     if ((transpose>2) || (transpose<0))
546     {     {
547      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");
548     }     }
549     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
550     {     {
551      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
552      Data ltemp(left);          Data ltemp(left);
553      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
554      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
555     }     }
556     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
557     {     {
558      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
559      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
560     }     }
561  //    left->operandCheck(*right);  //    left->operandCheck(*right);
562    
563     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
564     {     {
565      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
566     }     }
567     else     else
568     {     {
569      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
570     }     }
571     if (right->isLazy())     if (right->isLazy())
572     {     {
573      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
574     }     }
575     else     else
576     {     {
577      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
578     }     }
579     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
580     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
581     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
582     {     {
583      m_readytype='E';          m_readytype='E';
584     }     }
585     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
586     {     {
587      m_readytype='T';          m_readytype='T';
588     }     }
589     else     else
590     {     {
591      m_readytype='C';          m_readytype='C';
592     }     }
593     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);  
594     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
595     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  
596     LazyNodeSetup();     LazyNodeSetup();
 #endif  
597     SIZELIMIT     SIZELIMIT
598  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
599  }  }
600    
601    
602  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
603      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),          : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
604      m_op(op),          m_op(op),
605      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
606      m_transpose(0),          m_transpose(0),
607      m_tol(0)          m_tol(0)
608  {  {
609     if ((getOpgroup(op)!=G_NP1OUT_P))     if ((getOpgroup(op)!=G_NP1OUT_P))
610     {     {
611      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.");
612     }     }
613     DataLazy_ptr lleft;     DataLazy_ptr lleft;
614     if (!left->isLazy())     if (!left->isLazy())
615     {     {
616      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
617     }     }
618     else     else
619     {     {
620      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
621     }     }
622     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
623     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
624     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
625     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
626     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
627     LazyNodeSetup();     LazyNodeSetup();
 #endif  
628     SIZELIMIT     SIZELIMIT
629  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
630  }  }
631    
632  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
633      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
634      m_op(op),          m_op(op),
635      m_axis_offset(0),          m_axis_offset(0),
636      m_transpose(0),          m_transpose(0),
637      m_tol(tol)          m_tol(tol)
638  {  {
639     if ((getOpgroup(op)!=G_UNARY_P))     if ((getOpgroup(op)!=G_UNARY_P))
640     {     {
641      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.");
642     }     }
643     DataLazy_ptr lleft;     DataLazy_ptr lleft;
644     if (!left->isLazy())     if (!left->isLazy())
645     {     {
646      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
647     }     }
648     else     else
649     {     {
650      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
651     }     }
652     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
653     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
654     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
655     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
656     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
657     LazyNodeSetup();     LazyNodeSetup();
 #endif  
658     SIZELIMIT     SIZELIMIT
659  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
660  }  }
661    
662    
663  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)
664      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),          : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
665      m_op(op),          m_op(op),
666      m_axis_offset(axis0),          m_axis_offset(axis0),
667      m_transpose(axis1),          m_transpose(axis1),
668      m_tol(0)          m_tol(0)
669  {  {
670     if ((getOpgroup(op)!=G_NP1OUT_2P))     if ((getOpgroup(op)!=G_NP1OUT_2P))
671     {     {
672      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.");
673     }     }
674     DataLazy_ptr lleft;     DataLazy_ptr lleft;
675     if (!left->isLazy())     if (!left->isLazy())
676     {     {
677      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
678     }     }
679     else     else
680     {     {
681      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
682     }     }
683     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
684     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
685     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
686     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
687     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
 #ifdef LAZY_NODE_STORAGE  
688     LazyNodeSetup();     LazyNodeSetup();
 #endif  
689     SIZELIMIT     SIZELIMIT
690  LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
691  }  }
692    
693  DataLazy::~DataLazy()  
694    namespace
695  {  {
 #ifdef LAZY_NODE_SETUP  
    delete[] m_sampleids;  
    delete[] m_samples;  
 #endif  
 }  
696    
697        inline int max3(int a, int b, int c)
698        {
699            int t=(a>b?a:b);
700            return (t>c?t:c);
701    
702  int      }
 DataLazy::getBuffsRequired() const  
 {  
     return m_buffsRequired;  
703  }  }
704    
705    DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
706  size_t          : parent(left->getFunctionSpace(), left->getShape()),
707  DataLazy::getMaxSampleSize() const          m_op(CONDEVAL),
708            m_axis_offset(0),
709            m_transpose(0),
710            m_tol(0)
711  {  {
712      return m_maxsamplesize;  
713       DataLazy_ptr lmask;
714       DataLazy_ptr lleft;
715       DataLazy_ptr lright;
716       if (!mask->isLazy())
717       {
718            lmask=DataLazy_ptr(new DataLazy(mask));
719       }
720       else
721       {
722            lmask=dynamic_pointer_cast<DataLazy>(mask);
723       }
724       if (!left->isLazy())
725       {
726            lleft=DataLazy_ptr(new DataLazy(left));
727       }
728       else
729       {
730            lleft=dynamic_pointer_cast<DataLazy>(left);
731       }
732       if (!right->isLazy())
733       {
734            lright=DataLazy_ptr(new DataLazy(right));
735       }
736       else
737       {
738            lright=dynamic_pointer_cast<DataLazy>(right);
739       }
740       m_readytype=lmask->m_readytype;
741       if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
742       {
743            throw DataException("Programmer Error - condEval arguments must have the same readytype");
744       }
745       m_left=lleft;
746       m_right=lright;
747       m_mask=lmask;
748       m_samplesize=getNumDPPSample()*getNoValues();
749       m_children=m_left->m_children+m_right->m_children+m_mask->m_children+1;
750       m_height=max3(m_left->m_height,m_right->m_height,m_mask->m_height)+1;
751       LazyNodeSetup();
752       SIZELIMIT
753    LAZYDEBUG(cout << "(8)Lazy created with " << m_samplesize << endl;)
754  }  }
755    
756    
757    
758  size_t  DataLazy::~DataLazy()
 DataLazy::getSampleBufferSize() const  
759  {  {
760      return m_maxsamplesize*(max(1,m_buffsRequired));     delete[] m_sampleids;
761  }  }
762    
763    
764  /*  /*
765    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
766    This does the work for the collapse method.    This does the work for the collapse method.
767    For reasons of efficiency do not call this method on DataExpanded nodes.    For reasons of efficiency do not call this method on DataExpanded nodes.
768  */  */
769  DataReady_ptr  DataReady_ptr
770  DataLazy::collapseToReady()  DataLazy::collapseToReady() const
771  {  {
772    if (m_readytype=='E')    if (m_readytype=='E')
773    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
774      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
775    }    }
776    if (m_op==IDENTITY)    if (m_op==IDENTITY)
# Line 818  DataLazy::collapseToReady() Line 788  DataLazy::collapseToReady()
788    switch(m_op)    switch(m_op)
789    {    {
790      case ADD:      case ADD:
791      result=left+right;          result=left+right;
792      break;          break;
793      case SUB:            case SUB:          
794      result=left-right;          result=left-right;
795      break;          break;
796      case MUL:            case MUL:          
797      result=left*right;          result=left*right;
798      break;          break;
799      case DIV:            case DIV:          
800      result=left/right;          result=left/right;
801      break;          break;
802      case SIN:      case SIN:
803      result=left.sin();            result=left.sin();      
804      break;          break;
805      case COS:      case COS:
806      result=left.cos();          result=left.cos();
807      break;          break;
808      case TAN:      case TAN:
809      result=left.tan();          result=left.tan();
810      break;          break;
811      case ASIN:      case ASIN:
812      result=left.asin();          result=left.asin();
813      break;          break;
814      case ACOS:      case ACOS:
815      result=left.acos();          result=left.acos();
816      break;          break;
817      case ATAN:      case ATAN:
818      result=left.atan();          result=left.atan();
819      break;          break;
820      case SINH:      case SINH:
821      result=left.sinh();          result=left.sinh();
822      break;          break;
823      case COSH:      case COSH:
824      result=left.cosh();          result=left.cosh();
825      break;          break;
826      case TANH:      case TANH:
827      result=left.tanh();          result=left.tanh();
828      break;          break;
829      case ERF:      case ERF:
830      result=left.erf();          result=left.erf();
831      break;          break;
832     case ASINH:     case ASINH:
833      result=left.asinh();          result=left.asinh();
834      break;          break;
835     case ACOSH:     case ACOSH:
836      result=left.acosh();          result=left.acosh();
837      break;          break;
838     case ATANH:     case ATANH:
839      result=left.atanh();          result=left.atanh();
840      break;          break;
841      case LOG10:      case LOG10:
842      result=left.log10();          result=left.log10();
843      break;          break;
844      case LOG:      case LOG:
845      result=left.log();          result=left.log();
846      break;          break;
847      case SIGN:      case SIGN:
848      result=left.sign();          result=left.sign();
849      break;          break;
850      case ABS:      case ABS:
851      result=left.abs();          result=left.abs();
852      break;          break;
853      case NEG:      case NEG:
854      result=left.neg();          result=left.neg();
855      break;          break;
856      case POS:      case POS:
857      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
858      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
859      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
860      break;          break;
861      case EXP:      case EXP:
862      result=left.exp();          result=left.exp();
863      break;          break;
864      case SQRT:      case SQRT:
865      result=left.sqrt();          result=left.sqrt();
866      break;          break;
867      case RECIP:      case RECIP:
868      result=left.oneOver();          result=left.oneOver();
869      break;          break;
870      case GZ:      case GZ:
871      result=left.wherePositive();          result=left.wherePositive();
872      break;          break;
873      case LZ:      case LZ:
874      result=left.whereNegative();          result=left.whereNegative();
875      break;          break;
876      case GEZ:      case GEZ:
877      result=left.whereNonNegative();          result=left.whereNonNegative();
878      break;          break;
879      case LEZ:      case LEZ:
880      result=left.whereNonPositive();          result=left.whereNonPositive();
881      break;          break;
882      case NEZ:      case NEZ:
883      result=left.whereNonZero(m_tol);          result=left.whereNonZero(m_tol);
884      break;          break;
885      case EZ:      case EZ:
886      result=left.whereZero(m_tol);          result=left.whereZero(m_tol);
887      break;          break;
888      case SYM:      case SYM:
889      result=left.symmetric();          result=left.symmetric();
890      break;          break;
891      case NSYM:      case NSYM:
892      result=left.nonsymmetric();          result=left.antisymmetric();
893      break;          break;
894      case PROD:      case PROD:
895      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);          result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
896      break;          break;
897      case TRANS:      case TRANS:
898      result=left.transpose(m_axis_offset);          result=left.transpose(m_axis_offset);
899      break;          break;
900      case TRACE:      case TRACE:
901      result=left.trace(m_axis_offset);          result=left.trace(m_axis_offset);
902      break;          break;
903      case SWAP:      case SWAP:
904      result=left.swapaxes(m_axis_offset, m_transpose);          result=left.swapaxes(m_axis_offset, m_transpose);
905            break;
906        case MINVAL:
907            result=left.minval();
908            break;
909        case MAXVAL:
910            result=left.minval();
911            break;
912        case HER:
913        result=left.hermitian();
914      break;      break;
915      default:      default:
916      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)+".");
917    }    }
918    return result.borrowReadyPtr();    return result.borrowReadyPtr();
919  }  }
# Line 946  DataLazy::collapseToReady() Line 925  DataLazy::collapseToReady()
925     the purpose of using DataLazy in the first place).     the purpose of using DataLazy in the first place).
926  */  */
927  void  void
928  DataLazy::collapse()  DataLazy::collapse() const
929  {  {
930    if (m_op==IDENTITY)    if (m_op==IDENTITY)
931    {    {
932      return;          return;
933    }    }
934    if (m_readytype=='E')    if (m_readytype=='E')
935    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
936      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
937    }    }
938    m_id=collapseToReady();    m_id=collapseToReady();
939    m_op=IDENTITY;    m_op=IDENTITY;
940  }  }
941    
 /*  
   \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->resolveSample(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 (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->resolveSample(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->resolveSample(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->resolveSample(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->resolveSample(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->resolveSample(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->resolveSample(v,gap,sampleNo,lroffset);  
   gap+=m_left->getMaxSampleSize();  
   
   
 LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)  
   
   
   const ValueType* right=m_right->resolveSample(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  
   
942  // The result will be stored in m_samples  // The result will be stored in m_samples
943  // 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
944  const DataTypes::ValueType*  const DataTypes::RealVectorType*
945  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
946  {  {
947  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
948      // 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
949    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
950    {    {
951      collapse();          collapse();
952    }    }
953    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
954    {    {
955      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);  
 //     }  
956      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
957    #ifdef LAZY_STACK_PROF
958    int x;
959    if (&x<stackend[omp_get_thread_num()])
960    {
961           stackend[omp_get_thread_num()]=&x;
962    }
963    #endif
964      return &(vec);      return &(vec);
965    }    }
966    if (m_readytype!='E')    if (m_readytype!='E')
# Line 1615  LAZYDEBUG(cout << "Resolve sample " << t Line 969  LAZYDEBUG(cout << "Resolve sample " << t
969    }    }
970    if (m_sampleids[tid]==sampleNo)    if (m_sampleids[tid]==sampleNo)
971    {    {
972      roffset=tid*m_samplesize;          roffset=tid*m_samplesize;
973      return &(m_samples);        // sample is already resolved          return &(m_samples_r);            // sample is already resolved
974    }    }
975    m_sampleids[tid]=sampleNo;    m_sampleids[tid]=sampleNo;
976    
977    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
978    {    {
979    case G_UNARY:    case G_UNARY:
# Line 1628  LAZYDEBUG(cout << "Resolve sample " << t Line 983  LAZYDEBUG(cout << "Resolve sample " << t
983    case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);    case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
984    case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);    case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
985    case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);    case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
986      case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
987      case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
988    default:    default:
989      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)+".");
990    }    }
991  }  }
992    
993  const DataTypes::ValueType*  const DataTypes::RealVectorType*
994  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
995  {  {
996      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
997      // 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
998      // processing single points.          // processing single points.
999      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1000    if (m_readytype!='E')    if (m_readytype!='E')
1001    {    {
1002      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 1648  DataLazy::resolveNodeUnary(int tid, int Line 1005  DataLazy::resolveNodeUnary(int tid, int
1005    {    {
1006      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1007    }    }
1008    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1009    const double* left=&((*leftres)[roffset]);    const double* left=&((*leftres)[roffset]);
1010    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1011    double* result=&(m_samples[roffset]);    double* result=&(m_samples_r[roffset]);
1012    switch (m_op)    if (m_op==POS)
1013    {    {
1014      case SIN:        // this should be prevented earlier
1015      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);      // operation is meaningless for lazy
1016      break;          throw DataException("Programmer error - POS not supported for lazy data.");    
1017      case COS:    }
1018      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);    tensor_unary_array_operation(m_samplesize,
1019      break;                               left,
1020      case TAN:                               result,
1021      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);                               m_op,
1022      break;                               m_tol);  
1023      case ASIN:    return &(m_samples_r);
1024      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);  }
1025      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;  
1026    
1027    const DataTypes::RealVectorType*
1028    DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1029    {
1030            // we assume that any collapsing has been done before we get here
1031            // since we only have one argument we don't need to think about only
1032            // processing single points.
1033            // we will also know we won't get identity nodes
1034      if (m_readytype!='E')
1035      {
1036        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1037      }
1038      if (m_op==IDENTITY)
1039      {
1040        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1041      }
1042      size_t loffset=0;
1043      const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1044    
1045      roffset=m_samplesize*tid;
1046      unsigned int ndpps=getNumDPPSample();
1047      unsigned int psize=DataTypes::noValues(m_left->getShape());
1048      double* result=&(m_samples_r[roffset]);
1049      switch (m_op)
1050      {
1051        case MINVAL:
1052            {
1053              for (unsigned int z=0;z<ndpps;++z)
1054              {
1055                FMin op;
1056                *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1057                loffset+=psize;
1058                result++;
1059              }
1060            }
1061            break;
1062        case MAXVAL:
1063            {
1064              for (unsigned int z=0;z<ndpps;++z)
1065              {
1066              FMax op;
1067              *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1068              loffset+=psize;
1069              result++;
1070              }
1071            }
1072            break;
1073      default:      default:
1074      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1075    }    }
1076    return &(m_samples);    return &(m_samples_r);
1077  }  }
1078    
1079    const DataTypes::RealVectorType*
1080  const DataTypes::ValueType*  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
 DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)  
1081  {  {
1082      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1083      // 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
1084      // processing single points.          // processing single points.
1085    if (m_readytype!='E')    if (m_readytype!='E')
1086    {    {
1087      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 1780  DataLazy::resolveNodeNP1OUT(int tid, int Line 1091  DataLazy::resolveNodeNP1OUT(int tid, int
1091      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1092    }    }
1093    size_t subroffset;    size_t subroffset;
1094    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1095    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1096    size_t loop=0;    size_t loop=0;
1097    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
# Line 1789  DataLazy::resolveNodeNP1OUT(int tid, int Line 1100  DataLazy::resolveNodeNP1OUT(int tid, int
1100    switch (m_op)    switch (m_op)
1101    {    {
1102      case SYM:      case SYM:
1103      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1104      {          {
1105          DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              escript::symmetric(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(), offset);
1106          subroffset+=step;              subroffset+=step;
1107          offset+=step;              offset+=step;
1108      }          }
1109      break;          break;
1110      case NSYM:      case NSYM:
1111      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1112      {          {
1113          DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              escript::antisymmetric(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(), offset);
1114          subroffset+=step;              subroffset+=step;
1115          offset+=step;              offset+=step;
1116      }          }
1117      break;          break;
1118      default:      default:
1119      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1120    }    }
1121    return &m_samples;    return &m_samples_r;
1122  }  }
1123    
1124  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1125  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1126  {  {
1127      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1128      // 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
1129      // processing single points.          // processing single points.
1130    if (m_readytype!='E')    if (m_readytype!='E')
1131    {    {
1132      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 1826  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1137  DataLazy::resolveNodeNP1OUT_P(int tid, i
1137    }    }
1138    size_t subroffset;    size_t subroffset;
1139    size_t offset;    size_t offset;
1140    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1141    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1142    offset=roffset;    offset=roffset;
1143    size_t loop=0;    size_t loop=0;
# Line 1836  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1147  DataLazy::resolveNodeNP1OUT_P(int tid, i
1147    switch (m_op)    switch (m_op)
1148    {    {
1149      case TRACE:      case TRACE:
1150      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1151      {          {
1152              DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);              escript::trace(*leftres,m_left->getShape(),subroffset, m_samples_r ,getShape(),offset,m_axis_offset);
1153          subroffset+=instep;              subroffset+=instep;
1154          offset+=outstep;              offset+=outstep;
1155      }          }
1156      break;          break;
1157      case TRANS:      case TRANS:
1158      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1159      {          {
1160              DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);              escript::transpose(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(),offset,m_axis_offset);
1161          subroffset+=instep;              subroffset+=instep;
1162          offset+=outstep;              offset+=outstep;
1163      }          }
1164      break;          break;
1165      default:      default:
1166      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1167    }    }
1168    return &m_samples;    return &m_samples_r;
1169  }  }
1170    
1171    
1172  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1173  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1174  {  {
1175    if (m_readytype!='E')    if (m_readytype!='E')
1176    {    {
# Line 1871  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1182  DataLazy::resolveNodeNP1OUT_2P(int tid,
1182    }    }
1183    size_t subroffset;    size_t subroffset;
1184    size_t offset;    size_t offset;
1185    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1186    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1187    offset=roffset;    offset=roffset;
1188    size_t loop=0;    size_t loop=0;
# Line 1881  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1192  DataLazy::resolveNodeNP1OUT_2P(int tid,
1192    switch (m_op)    switch (m_op)
1193    {    {
1194      case SWAP:      case SWAP:
1195      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1196      {          {
1197              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_r, getShape(),offset, m_axis_offset, m_transpose);
1198          subroffset+=instep;              subroffset+=instep;
1199          offset+=outstep;              offset+=outstep;
1200      }          }
1201      break;          break;
1202      default:      default:
1203      throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1204    }    }
1205    return &m_samples;    return &m_samples_r;
1206  }  }
1207    
1208    const DataTypes::RealVectorType*
1209    DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1210    {
1211      if (m_readytype!='E')
1212      {
1213        throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1214      }
1215      if (m_op!=CONDEVAL)
1216      {
1217        throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1218      }
1219      size_t subroffset;
1220    
1221      const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1222      const RealVectorType* srcres=0;
1223      if ((*maskres)[subroffset]>0)
1224      {
1225            srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1226      }
1227      else
1228      {
1229            srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1230      }
1231    
1232      // Now we need to copy the result
1233    
1234      roffset=m_samplesize*tid;
1235      for (int i=0;i<m_samplesize;++i)
1236      {
1237            m_samples_r[roffset+i]=(*srcres)[subroffset+i];  
1238      }
1239    
1240      return &m_samples_r;
1241    }
1242    
1243  // 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
1244  // 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 1905  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1249  DataLazy::resolveNodeNP1OUT_2P(int tid,
1249  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1250  // For example, 2+Vector.  // For example, 2+Vector.
1251  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1252  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1253  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1254  {  {
1255  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1256    
1257    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
1258      // first work out which of the children are expanded          // first work out which of the children are expanded
1259    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1260    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1261    if (!leftExp && !rightExp)    if (!leftExp && !rightExp)
1262    {    {
1263      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'.");
1264    }    }
1265    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1266    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1267    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1268    {    {
1269      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.");
1270    }    }
1271    size_t leftsize=m_left->getNoValues();    size_t leftsize=m_left->getNoValues();
1272    size_t rightsize=m_right->getNoValues();    size_t rightsize=m_right->getNoValues();
1273    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
1274    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
1275    int rightstep=0;    int rightstep=0;
1276    int numsteps=0;       // total number of steps for the inner loop    int numsteps=0;               // total number of steps for the inner loop
1277    int oleftstep=0;  // the o variables refer to the outer loop    int oleftstep=0;      // the o variables refer to the outer loop
1278    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
1279    int onumsteps=1;    int onumsteps=1;
1280        
1281    bool LES=(leftExp && leftScalar); // Left is an expanded scalar    bool LES=(leftExp && leftScalar);     // Left is an expanded scalar
1282    bool RES=(rightExp && rightScalar);    bool RES=(rightExp && rightScalar);
1283    bool LS=(!leftExp && leftScalar); // left is a single scalar    bool LS=(!leftExp && leftScalar);     // left is a single scalar
1284    bool RS=(!rightExp && rightScalar);    bool RS=(!rightExp && rightScalar);
1285    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1286    bool RN=(!rightExp && !rightScalar);    bool RN=(!rightExp && !rightScalar);
1287    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1288    bool REN=(rightExp && !rightScalar);    bool REN=(rightExp && !rightScalar);
1289    
1290    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
1291    {    {
1292      chunksize=m_left->getNumDPPSample()*leftsize;          chunksize=m_left->getNumDPPSample()*leftsize;
1293      leftstep=0;          leftstep=0;
1294      rightstep=0;          rightstep=0;
1295      numsteps=1;          numsteps=1;
1296    }    }
1297    else if (LES || RES)    else if (LES || RES)
1298    {    {
1299      chunksize=1;          chunksize=1;
1300      if (LES)        // left is an expanded scalar          if (LES)                // left is an expanded scalar
1301      {          {
1302          if (RS)                  if (RS)
1303          {                  {
1304             leftstep=1;                     leftstep=1;
1305             rightstep=0;                     rightstep=0;
1306             numsteps=m_left->getNumDPPSample();                     numsteps=m_left->getNumDPPSample();
1307          }                  }
1308          else        // RN or REN                  else            // RN or REN
1309          {                  {
1310             leftstep=0;                     leftstep=0;
1311             oleftstep=1;                     oleftstep=1;
1312             rightstep=1;                     rightstep=1;
1313             orightstep=(RN ? -(int)rightsize : 0);                     orightstep=(RN ? -(int)rightsize : 0);
1314             numsteps=rightsize;                     numsteps=rightsize;
1315             onumsteps=m_left->getNumDPPSample();                     onumsteps=m_left->getNumDPPSample();
1316          }                  }
1317      }          }
1318      else        // right is an expanded scalar          else            // right is an expanded scalar
1319      {          {
1320          if (LS)                  if (LS)
1321          {                  {
1322             rightstep=1;                     rightstep=1;
1323             leftstep=0;                     leftstep=0;
1324             numsteps=m_right->getNumDPPSample();                     numsteps=m_right->getNumDPPSample();
1325          }                  }
1326          else                  else
1327          {                  {
1328             rightstep=0;                     rightstep=0;
1329             orightstep=1;                     orightstep=1;
1330             leftstep=1;                     leftstep=1;
1331             oleftstep=(LN ? -(int)leftsize : 0);                     oleftstep=(LN ? -(int)leftsize : 0);
1332             numsteps=leftsize;                     numsteps=leftsize;
1333             onumsteps=m_right->getNumDPPSample();                     onumsteps=m_right->getNumDPPSample();
1334          }                  }
1335      }          }
1336    }    }
1337    else  // this leaves (LEN, RS), (LEN, RN) and their transposes    else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1338    {    {
1339      if (LEN)    // and Right will be a single value          if (LEN)        // and Right will be a single value
1340      {          {
1341          chunksize=rightsize;                  chunksize=rightsize;
1342          leftstep=rightsize;                  leftstep=rightsize;
1343          rightstep=0;                  rightstep=0;
1344          numsteps=m_left->getNumDPPSample();                  numsteps=m_left->getNumDPPSample();
1345          if (RS)                  if (RS)
1346          {                  {
1347             numsteps*=leftsize;                     numsteps*=leftsize;
1348          }                  }
1349      }          }
1350      else    // REN          else    // REN
1351      {          {
1352          chunksize=leftsize;                  chunksize=leftsize;
1353          rightstep=leftsize;                  rightstep=leftsize;
1354          leftstep=0;                  leftstep=0;
1355          numsteps=m_right->getNumDPPSample();                  numsteps=m_right->getNumDPPSample();
1356          if (LS)                  if (LS)
1357          {                  {
1358             numsteps*=rightsize;                     numsteps*=rightsize;
1359          }                  }
1360      }          }
1361    }    }
1362    
1363    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1364      // Get the values of sub-expressions          // Get the values of sub-expressions
1365    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1366    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1367  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1368  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;)
1369  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
# Line 2033  LAZYDEBUG(cout << "Right res["<< rroffse Line 1377  LAZYDEBUG(cout << "Right res["<< rroffse
1377    
1378    
1379    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1380    double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved    double* resultp=&(m_samples_r[roffset]);                // results are stored at the vector offset we received
1381    switch(m_op)    switch(m_op)
1382    {    {
1383      case ADD:      case ADD:
1384          PROC_OP(NO_ARG,plus<double>());          //PROC_OP(NO_ARG,plus<double>());
1385      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1386                 &(*left)[0],
1387                 &(*right)[0],
1388                 chunksize,
1389                 onumsteps,
1390                 numsteps,
1391                 resultStep,
1392                 leftstep,
1393                 rightstep,
1394                 oleftstep,
1395                 orightstep,
1396                 lroffset,
1397                 rroffset,
1398                 escript::ES_optype::ADD);  
1399            break;
1400      case SUB:      case SUB:
1401      PROC_OP(NO_ARG,minus<double>());        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1402      break;               &(*left)[0],
1403                 &(*right)[0],
1404                 chunksize,
1405                 onumsteps,
1406                 numsteps,
1407                 resultStep,
1408                 leftstep,
1409                 rightstep,
1410                 oleftstep,
1411                 orightstep,
1412                 lroffset,
1413                 rroffset,
1414                 escript::ES_optype::SUB);        
1415            //PROC_OP(NO_ARG,minus<double>());
1416            break;
1417      case MUL:      case MUL:
1418      PROC_OP(NO_ARG,multiplies<double>());          //PROC_OP(NO_ARG,multiplies<double>());
1419      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1420                 &(*left)[0],
1421                 &(*right)[0],
1422                 chunksize,
1423                 onumsteps,
1424                 numsteps,
1425                 resultStep,
1426                 leftstep,
1427                 rightstep,
1428                 oleftstep,
1429                 orightstep,
1430                 lroffset,
1431                 rroffset,
1432                 escript::ES_optype::MUL);        
1433            break;
1434      case DIV:      case DIV:
1435      PROC_OP(NO_ARG,divides<double>());          //PROC_OP(NO_ARG,divides<double>());
1436      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1437                 &(*left)[0],
1438                 &(*right)[0],
1439                 chunksize,
1440                 onumsteps,
1441                 numsteps,
1442                 resultStep,
1443                 leftstep,
1444                 rightstep,
1445                 oleftstep,
1446                 orightstep,
1447                 lroffset,
1448                 rroffset,
1449                 escript::ES_optype::DIV);        
1450            break;
1451      case POW:      case POW:
1452         PROC_OP(double (double,double),::pow);         //PROC_OP(double (double,double),::pow);
1453      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1454                 &(*left)[0],
1455                 &(*right)[0],
1456                 chunksize,
1457                 onumsteps,
1458                 numsteps,
1459                 resultStep,
1460                 leftstep,
1461                 rightstep,
1462                 oleftstep,
1463                 orightstep,
1464                 lroffset,
1465                 rroffset,
1466                 escript::ES_optype::POW);        
1467            break;
1468      default:      default:
1469      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1470    }    }
1471  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples_r[roffset] << endl;)
1472    return &m_samples;    return &m_samples_r;
1473  }  }
1474    
1475    
1476  // 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
1477  // 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.
1478  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1479  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1480  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1481  {  {
1482  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1483    
1484    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
1485      // first work out which of the children are expanded          // first work out which of the children are expanded
1486    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1487    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1488    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1489    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
1490    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1491    
1492    int resultStep=getNoValues();    int resultStep=getNoValues();
1493    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1494    size_t offset=roffset;    size_t offset=roffset;
1495    
1496    const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);    const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1497    
1498    const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1499    
1500  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;
1501  cout << getNoValues() << endl;)  cout << getNoValues() << endl;)
# Line 2095  LAZYDEBUG(cout << "m_samplesize=" << m_s Line 1509  LAZYDEBUG(cout << "m_samplesize=" << m_s
1509  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1510  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1511    
1512    double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(m_samples_r[offset]);         // results are stored at the vector offset we received
1513    switch(m_op)    switch(m_op)
1514    {    {
1515      case PROD:      case PROD:
1516      for (int i=0;i<steps;++i,resultp+=resultStep)          for (int i=0;i<steps;++i,resultp+=resultStep)
1517      {          {
1518            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1519            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1520    
1521  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1522  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1523    
1524            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);
1525    
1526        lroffset+=leftStep;            lroffset+=leftStep;
1527        rroffset+=rightStep;            rroffset+=rightStep;
1528      }          }
1529      break;          break;
1530      default:      default:
1531      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1532    }    }
1533    roffset=offset;    roffset=offset;
1534    return &m_samples;    return &m_samples_r;
1535  }  }
 #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::resolveSample(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);  
   default:  
     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
   }  
1536    
 }  
1537    
1538  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1539  DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)  DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1540  {  {
1541  #ifdef _OPENMP  #ifdef _OPENMP
1542      int tid=omp_get_thread_num();          int tid=omp_get_thread_num();
1543  #else  #else
1544      int tid=0;          int tid=0;
1545  #endif  #endif
1546      return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);  
1547    #ifdef LAZY_STACK_PROF
1548            stackstart[tid]=&tid;
1549            stackend[tid]=&tid;
1550            const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1551            size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1552            #pragma omp critical
1553            if (d>maxstackuse)
1554            {
1555    cout << "Max resolve Stack use " << d << endl;
1556                    maxstackuse=d;
1557            }
1558            return r;
1559    #else
1560            return resolveNodeSample(tid, sampleNo, roffset);
1561    #endif
1562  }  }
1563    
1564    
# Line 2195  void Line 1567  void
1567  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
1568  {  {
1569     if (m_op==IDENTITY)     if (m_op==IDENTITY)
1570      return;          return;
 #ifndef LAZY_NODE_STORAGE  
    DataReady_ptr p=resolveVectorWorker();  
 #else  
1571     DataReady_ptr p=resolveNodeWorker();     DataReady_ptr p=resolveNodeWorker();
 #endif  
1572     makeIdentity(p);     makeIdentity(p);
1573  }  }
1574    
# Line 2216  void DataLazy::makeIdentity(const DataRe Line 1584  void DataLazy::makeIdentity(const DataRe
1584     else if(p->isExpanded()) {m_readytype='E';}     else if(p->isExpanded()) {m_readytype='E';}
1585     else if (p->isTagged()) {m_readytype='T';}     else if (p->isTagged()) {m_readytype='T';}
1586     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
    m_buffsRequired=1;  
1587     m_samplesize=p->getNumDPPSample()*p->getNoValues();     m_samplesize=p->getNumDPPSample()*p->getNoValues();
    m_maxsamplesize=m_samplesize;  
1588     m_left.reset();     m_left.reset();
1589     m_right.reset();     m_right.reset();
1590       m_iscompl=p->isComplex();
1591  }  }
1592    
1593    
# Line 2231  DataLazy::resolve() Line 1598  DataLazy::resolve()
1598      return m_id;      return m_id;
1599  }  }
1600    
 #ifdef LAZY_NODE_STORAGE  
1601    
1602  // 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 */
1603  DataReady_ptr  void
1604  DataLazy::resolveNodeWorker()  DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
1605  {  {
1606    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally    if (dats.empty())
1607    {    {
1608      collapse();          return;
1609    }    }
1610    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.    vector<DataLazy*> work;
1611    {    FunctionSpace fs=dats[0]->getFunctionSpace();
1612      return m_id;    bool match=true;
1613    }    for (int i=dats.size()-1;i>=0;--i)
1614      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'    {
1615    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));          if (dats[i]->m_readytype!='E')
1616    ValueType& resvec=result->getVectorRW();          {
1617    DataReady_ptr resptr=DataReady_ptr(result);                  dats[i]->collapse();
1618            }
1619    int sample;          if (dats[i]->m_op!=IDENTITY)
1620    int totalsamples=getNumSamples();          {
1621    const ValueType* res=0;   // Storage for answer                  work.push_back(dats[i]);
1622  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)                  if (fs!=dats[i]->getFunctionSpace())
1623    #pragma omp parallel for private(sample,res) schedule(static)                  {
1624    for (sample=0;sample<totalsamples;++sample)                          match=false;
1625    {                  }
1626      size_t roffset=0;          }
1627      }
1628      if (work.empty())
1629      {
1630            return;         // no work to do
1631      }
1632      if (match)    // all functionspaces match.  Yes I realise this is overly strict
1633      {             // it is possible that dats[0] is one of the objects which we discarded and
1634                    // all the other functionspaces match.
1635            vector<DataExpanded*> dep;
1636            vector<RealVectorType*> vecs;
1637            for (int i=0;i<work.size();++i)
1638            {
1639                    dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
1640                    vecs.push_back(&(dep[i]->getVectorRW()));
1641            }
1642            int totalsamples=work[0]->getNumSamples();
1643            const RealVectorType* res=0; // Storage for answer
1644            int sample;
1645            #pragma omp parallel private(sample, res)
1646            {
1647                size_t roffset=0;
1648                #pragma omp for schedule(static)
1649                for (sample=0;sample<totalsamples;++sample)
1650                {
1651                    roffset=0;
1652                    int j;
1653                    for (j=work.size()-1;j>=0;--j)
1654                    {
1655  #ifdef _OPENMP  #ifdef _OPENMP
1656      res=resolveNodeSample(omp_get_thread_num(),sample,roffset);                      res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1657  #else  #else
1658      res=resolveNodeSample(0,sample,roffset);                      res=work[j]->resolveNodeSample(0,sample,roffset);
1659  #endif  #endif
1660  LAZYDEBUG(cout << "Sample #" << sample << endl;)                      RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1661  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )                      memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
1662      DataVector::size_type outoffset=result->getPointOffset(sample,0);                  }
1663      memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));              }
1664            }
1665            // Now we need to load the new results as identity ops into the lazy nodes
1666            for (int i=work.size()-1;i>=0;--i)
1667            {
1668                work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1669            }
1670      }
1671      else  // functionspaces do not match
1672      {
1673            for (int i=0;i<work.size();++i)
1674            {
1675                    work[i]->resolveToIdentity();
1676            }
1677    }    }
   return resptr;  
1678  }  }
1679    
 #endif // LAZY_NODE_STORAGE  
1680    
1681  // To simplify the memory management, all threads operate on one large vector, rather than one each.  
1682  // Each sample is evaluated independently and copied into the result DataExpanded.  // This version of resolve uses storage in each node to hold results
1683  DataReady_ptr  DataReady_ptr
1684  DataLazy::resolveVectorWorker()  DataLazy::resolveNodeWorker()
1685  {  {
1686      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  
1687    {    {
1688      collapse();      collapse();
1689    }    }
1690    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.
1691    {    {
1692      return m_id;      return m_id;
1693    }    }
1694      // 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'
1695    size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  RealVectorType(getNoValues()));
1696      // 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();  
1697    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1698    
1699    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1700    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1701    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  
1702  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1703    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel private(sample,res)
   for (sample=0;sample<totalsamples;++sample)  
1704    {    {
1705  LAZYDEBUG(cout << "################################# " << sample << endl;)          size_t roffset=0;
1706    #ifdef LAZY_STACK_PROF
1707            stackstart[omp_get_thread_num()]=&roffset;
1708            stackend[omp_get_thread_num()]=&roffset;
1709    #endif
1710            #pragma omp for schedule(static)
1711            for (sample=0;sample<totalsamples;++sample)
1712            {
1713                    roffset=0;
1714  #ifdef _OPENMP  #ifdef _OPENMP
1715      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);                  res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1716  #else  #else
1717      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.                  res=resolveNodeSample(0,sample,roffset);
1718  #endif  #endif
1719  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1720  LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1721      outoffset=result->getPointOffset(sample,0);                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1722  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)                  memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1723      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector          }
1724      {    }
1725  LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)  #ifdef LAZY_STACK_PROF
1726      resvec[outoffset]=(*res)[resoffset];    for (int i=0;i<getNumberOfThreads();++i)
1727      }    {
1728  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]);
1729  LAZYDEBUG(cerr << "*********************************" << endl;)  //      cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1730            if (r>maxstackuse)
1731            {
1732                    maxstackuse=r;
1733            }
1734    }    }
1735      cout << "Max resolve Stack use=" << maxstackuse << endl;
1736    #endif
1737    return resptr;    return resptr;
1738  }  }
1739    
# Line 2335  std::string Line 1741  std::string
1741  DataLazy::toString() const  DataLazy::toString() const
1742  {  {
1743    ostringstream oss;    ostringstream oss;
1744    oss << "Lazy Data:";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1745    intoString(oss);    switch (escriptParams.getLazyStrFmt())
1746      {
1747      case 1:       // tree format
1748            oss << endl;
1749            intoTreeString(oss,"");
1750            break;
1751      case 2:       // just the depth
1752            break;
1753      default:
1754            intoString(oss);
1755            break;
1756      }
1757    return oss.str();    return oss.str();
1758  }  }
1759    
# Line 2348  DataLazy::intoString(ostringstream& oss) Line 1765  DataLazy::intoString(ostringstream& oss)
1765    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1766    {    {
1767    case G_IDENTITY:    case G_IDENTITY:
1768      if (m_id->isExpanded())          if (m_id->isExpanded())
1769      {          {
1770         oss << "E";             oss << "E";
1771      }          }
1772      else if (m_id->isTagged())          else if (m_id->isTagged())
1773      {          {
1774        oss << "T";            oss << "T";
1775      }          }
1776      else if (m_id->isConstant())          else if (m_id->isConstant())
1777      {          {
1778        oss << "C";            oss << "C";
1779      }          }
1780      else          else
1781      {          {
1782        oss << "?";            oss << "?";
1783      }          }
1784      oss << '@' << m_id.get();          oss << '@' << m_id.get();
1785      break;          break;
1786    case G_BINARY:    case G_BINARY:
1787      oss << '(';          oss << '(';
1788      m_left->intoString(oss);          m_left->intoString(oss);
1789      oss << ' ' << opToString(m_op) << ' ';          oss << ' ' << opToString(m_op) << ' ';
1790      m_right->intoString(oss);          m_right->intoString(oss);
1791      oss << ')';          oss << ')';
1792      break;          break;
1793    case G_UNARY:    case G_UNARY:
1794    case G_UNARY_P:    case G_UNARY_P:
1795    case G_NP1OUT:    case G_NP1OUT:
1796    case G_NP1OUT_P:    case G_NP1OUT_P:
1797      oss << opToString(m_op) << '(';    case G_REDUCTION:
1798      m_left->intoString(oss);          oss << opToString(m_op) << '(';
1799      oss << ')';          m_left->intoString(oss);
1800      break;          oss << ')';
1801            break;
1802    case G_TENSORPROD:    case G_TENSORPROD:
1803      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1804      m_left->intoString(oss);          m_left->intoString(oss);
1805      oss << ", ";          oss << ", ";
1806      m_right->intoString(oss);          m_right->intoString(oss);
1807      oss << ')';          oss << ')';
1808      break;          break;
1809    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1810      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1811      m_left->intoString(oss);          m_left->intoString(oss);
1812      oss << ", " << m_axis_offset << ", " << m_transpose;          oss << ", " << m_axis_offset << ", " << m_transpose;
1813      oss << ')';          oss << ')';
1814      break;          break;
1815      case G_CONDEVAL:
1816            oss << opToString(m_op)<< '(' ;
1817            m_mask->intoString(oss);
1818            oss << " ? ";
1819            m_left->intoString(oss);
1820            oss << " : ";
1821            m_right->intoString(oss);
1822            oss << ')';
1823            break;
1824    default:    default:
1825      oss << "UNKNOWN";          oss << "UNKNOWN";
1826    }    }
1827  }  }
1828    
1829    
1830    void
1831    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1832    {
1833      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1834      switch (getOpgroup(m_op))
1835      {
1836      case G_IDENTITY:
1837            if (m_id->isExpanded())
1838            {
1839               oss << "E";
1840            }
1841            else if (m_id->isTagged())
1842            {
1843              oss << "T";
1844            }
1845            else if (m_id->isConstant())
1846            {
1847              oss << "C";
1848            }
1849            else
1850            {
1851              oss << "?";
1852            }
1853            oss << '@' << m_id.get() << endl;
1854            break;
1855      case G_BINARY:
1856            oss << opToString(m_op) << endl;
1857            indent+='.';
1858            m_left->intoTreeString(oss, indent);
1859            m_right->intoTreeString(oss, indent);
1860            break;
1861      case G_UNARY:
1862      case G_UNARY_P:
1863      case G_NP1OUT:
1864      case G_NP1OUT_P:
1865      case G_REDUCTION:
1866            oss << opToString(m_op) << endl;
1867            indent+='.';
1868            m_left->intoTreeString(oss, indent);
1869            break;
1870      case G_TENSORPROD:
1871            oss << opToString(m_op) << endl;
1872            indent+='.';
1873            m_left->intoTreeString(oss, indent);
1874            m_right->intoTreeString(oss, indent);
1875            break;
1876      case G_NP1OUT_2P:
1877            oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1878            indent+='.';
1879            m_left->intoTreeString(oss, indent);
1880            break;
1881      default:
1882            oss << "UNKNOWN";
1883      }
1884    }
1885    
1886    
1887  DataAbstract*  DataAbstract*
1888  DataLazy::deepCopy()  DataLazy::deepCopy() const
1889  {  {
1890    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1891    {    {
1892    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1893    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY:
1894    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1895      case G_UNARY_P:       return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1896      case G_BINARY:        return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1897    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);
1898    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);
1899      case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
1900      case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1901    default:    default:
1902      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)+".");
1903    }    }
1904  }  }
1905    
1906    // For this, we don't care what op you were doing because the answer is now zero
1907    DataAbstract*
1908    DataLazy::zeroedCopy() const
1909    {
1910      return new DataLazy(m_id->zeroedCopy()->getPtr());
1911    }
1912    
1913  // There is no single, natural interpretation of getLength on DataLazy.  // There is no single, natural interpretation of getLength on DataLazy.
1914  // Instances of DataReady can look at the size of their vectors.  // Instances of DataReady can look at the size of their vectors.
# Line 2421  DataLazy::deepCopy() Line 1916  DataLazy::deepCopy()
1916  // 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
1917  // form part of the expression.  // form part of the expression.
1918  // Rather than have people making assumptions, I have disabled the method.  // Rather than have people making assumptions, I have disabled the method.
1919  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1920  DataLazy::getLength() const  DataLazy::getLength() const
1921  {  {
1922    throw DataException("getLength() does not make sense for lazy data.");    throw DataException("getLength() does not make sense for lazy data.");
# Line 2436  DataLazy::getSlice(const DataTypes::Regi Line 1931  DataLazy::getSlice(const DataTypes::Regi
1931    
1932    
1933  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
1934  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1935  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1936                   int dataPointNo)                   int dataPointNo)
1937  {  {
1938    if (m_op==IDENTITY)    if (m_op==IDENTITY)
1939    {    {
1940      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1941    }    }
1942    if (m_readytype!='E')    if (m_readytype!='E')
1943    {    {
1944      collapse();          collapse();
1945      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1946    }    }
1947    // 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
1948    // so we only need to know which child to ask    // so we only need to know which child to ask
1949    if (m_left->m_readytype=='E')    if (m_left->m_readytype=='E')
1950    {    {
1951      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
1952    }    }
1953    else    else
1954    {    {
1955      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
1956    }    }
1957  }  }
1958    
1959  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
1960  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1961  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1962                   int dataPointNo) const                   int dataPointNo) const
1963  {  {
1964    if (m_op==IDENTITY)    if (m_op==IDENTITY)
1965    {    {
1966      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1967    }    }
1968    if (m_readytype=='E')    if (m_readytype=='E')
1969    {    {
# Line 2476  DataLazy::getPointOffset(int sampleNo, Line 1971  DataLazy::getPointOffset(int sampleNo,
1971      // so we only need to know which child to ask      // so we only need to know which child to ask
1972      if (m_left->m_readytype=='E')      if (m_left->m_readytype=='E')
1973      {      {
1974      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
1975      }      }
1976      else      else
1977      {      {
1978      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
1979      }      }
1980    }    }
1981    if (m_readytype=='C')    if (m_readytype=='C')
1982    {    {
1983      return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter          return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1984    }    }
1985    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).");
1986  }  }
# Line 2495  DataLazy::getPointOffset(int sampleNo, Line 1990  DataLazy::getPointOffset(int sampleNo,
1990  void  void
1991  DataLazy::setToZero()  DataLazy::setToZero()
1992  {  {
1993  //   DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::RealVectorType v(getNoValues(),0);
1994  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1995  //   m_op=IDENTITY;  //   m_op=IDENTITY;
1996  //   m_right.reset();    //   m_right.reset();  
# Line 2503  DataLazy::setToZero() Line 1998  DataLazy::setToZero()
1998  //   m_readytype='C';  //   m_readytype='C';
1999  //   m_buffsRequired=1;  //   m_buffsRequired=1;
2000    
2001    privdebug=privdebug;  // to stop the compiler complaining about unused privdebug    (void)privdebug;  // to stop the compiler complaining about unused privdebug
2002    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).");
2003  }  }
2004    
2005  bool  bool
2006  DataLazy::actsExpanded() const  DataLazy::actsExpanded() const
2007  {  {
2008      return (m_readytype=='E');          return (m_readytype=='E');
2009  }  }
2010    
2011  }   // end namespace  } // end namespace
2012    

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

  ViewVC Help
Powered by ViewVC 1.1.26