/[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 4286 by caltinay, Thu Mar 7 04:28:11 2013 UTC trunk/escriptcore/src/DataLazy.cpp revision 6469 by jfenwick, Tue Jan 17 07:45:36 2017 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2016 by The University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
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)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * 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"
 #include "esysUtils/Esys_MPI.h"  
 #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  #ifdef USE_NETCDF  using namespace escript::DataTypes;
 #include <netcdfcpp.h>  
 #endif  
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 43  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())) {cerr << "SIZE LIMIT EXCEEDED " << m_height << endl;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 if (m_height>escript::escriptParams.getTOO_MANY_LEVELS())  {if (escript::escriptParams.getLAZY_VERBOSE()){cerr << "SIZE LIMIT EXCEEDED height=" << m_height << endl;}resolveToIdentity();}  #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?
# Line 69  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 large 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 132  std::vector<void*> stackend(getNumberOfT Line 133  std::vector<void*> stackend(getNumberOfT
133  size_t maxstackuse=0;  size_t maxstackuse=0;
134  #endif  #endif
135    
 enum ES_opgroup  
 {  
    G_UNKNOWN,  
    G_IDENTITY,  
    G_BINARY,        // pointwise operations with two arguments  
    G_UNARY,     // pointwise operations with one argument  
    G_UNARY_P,       // pointwise operations with one argument, requiring a parameter  
    G_NP1OUT,        // non-pointwise op with one output  
    G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter  
    G_TENSORPROD,    // general tensor product  
    G_NP1OUT_2P,     // non-pointwise op with one output requiring two params  
    G_REDUCTION,     // non-pointwise unary op with a scalar output  
    G_CONDEVAL  
 };  
   
   
   
   
 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",  
             "sin","cos","tan",  
             "asin","acos","atan","sinh","cosh","tanh","erf",  
             "asinh","acosh","atanh",  
             "log10","log","sign","abs","neg","pos","exp","sqrt",  
             "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",  
             "symmetric","nonsymmetric",  
             "prod",  
             "transpose", "trace",  
             "swapaxes",  
             "minval", "maxval",  
             "condEval"};  
 int ES_opcount=44;  
 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,  
             G_UNARY,G_UNARY,G_UNARY, //10  
             G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17  
             G_UNARY,G_UNARY,G_UNARY,                    // 20  
             G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 28  
             G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35  
             G_NP1OUT,G_NP1OUT,  
             G_TENSORPROD,  
             G_NP1OUT_P, G_NP1OUT_P,  
             G_NP1OUT_2P,  
             G_REDUCTION, G_REDUCTION,  
             G_CONDEVAL};  
 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();
# Line 198  resultFS(DataAbstract_ptr left, DataAbst Line 150  resultFS(DataAbstract_ptr left, DataAbst
150      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
151      if (res==1)      if (res==1)
152      {      {
153      return l;          return l;
154      }      }
155      if (res==-1)      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 214  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    
176        if (left->getRank()==0)   // we need to allow scalar * anything            if (left->getRank()==0)       // we need to allow scalar * anything
177        {            {
178          return right->getShape();                  return right->getShape();
179        }            }
180        if (right->getRank()==0)            if (right->getRank()==0)
181        {            {
182          return left->getShape();                  return left->getShape();
183        }            }
184        throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");            throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
185      }          }
186      return left->getShape();          return left->getShape();
187  }  }
188    
189  // return the shape for "op left"  // return the shape for "op left"
# Line 239  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              stringstream e;              stringstream e;
204              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
205              throw DataException(e.str());              throw DataException(e.str());
206          }          }
207          for (int i=0; i<rank; i++)          for (int i=0; i<rank; i++)
208          {                  {
209             int index = (axis_offset+i)%rank;                     int index = (axis_offset+i)%rank;
210             sh.push_back(s[index]); // Append to new shape             sh.push_back(s[index]); // Append to new shape
211          }          }
212          return sh;                  return sh;
213         }             }
214      break;          break;
215      case TRACE:          case TRACE:
216         {             {
217          int rank=left->getRank();                  int rank=left->getRank();
218          if (rank<2)                  if (rank<2)
219          {                  {
220             throw DataException("Trace can only be computed for objects with rank 2 or greater.");                     throw DataException("Trace can only be computed for objects with rank 2 or greater.");
221          }                  }
222          if ((axis_offset>rank-2) || (axis_offset<0))                  if ((axis_offset>rank-2) || (axis_offset<0))
223          {                  {
224             throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");                     throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
225          }                  }
226          if (rank==2)                  if (rank==2)
227          {                  {
228             return DataTypes::scalarShape;                     return DataTypes::scalarShape;
229          }                  }
230          else if (rank==3)                  else if (rank==3)
231          {                  {
232             DataTypes::ShapeType sh;                     DataTypes::ShapeType sh;
233                 if (axis_offset==0)                     if (axis_offset==0)
234             {                     {
235                  sh.push_back(left->getShape()[2]);                          sh.push_back(left->getShape()[2]);
236                 }                     }
237                 else     // offset==1                     else         // offset==1
238             {                     {
239              sh.push_back(left->getShape()[0]);                          sh.push_back(left->getShape()[0]);
240                 }                     }
241             return sh;                     return sh;
242          }                  }
243          else if (rank==4)                  else if (rank==4)
244          {                  {
245             DataTypes::ShapeType sh;                     DataTypes::ShapeType sh;
246             const DataTypes::ShapeType& s=left->getShape();                     const DataTypes::ShapeType& s=left->getShape();
247                 if (axis_offset==0)                     if (axis_offset==0)
248             {                     {
249                  sh.push_back(s[2]);                          sh.push_back(s[2]);
250                  sh.push_back(s[3]);                          sh.push_back(s[3]);
251                 }                     }
252                 else if (axis_offset==1)                     else if (axis_offset==1)
253             {                     {
254                  sh.push_back(s[0]);                          sh.push_back(s[0]);
255                  sh.push_back(s[3]);                          sh.push_back(s[3]);
256                 }                     }
257             else     // offset==2                     else         // offset==2
258             {                     {
259              sh.push_back(s[0]);                          sh.push_back(s[0]);
260              sh.push_back(s[1]);                          sh.push_back(s[1]);
261             }                     }
262             return sh;                     return sh;
263          }                  }
264          else        // unknown rank                  else            // unknown rank
265          {                  {
266             throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");                     throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
267          }                  }
268         }             }
269      break;          break;
270          default:          default:
271      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");          throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
272      }          }
273  }  }
274    
275  DataTypes::ShapeType  DataTypes::ShapeType
# Line 373  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 382  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 431  GTPShape(DataAbstract_ptr left, DataAbst Line 383  GTPShape(DataAbstract_ptr left, DataAbst
383    return shape2;    return shape2;
384  }  }
385    
386  }   // end anonymous namespace  }       // 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    
388  void DataLazy::LazyNodeSetup()  void DataLazy::LazyNodeSetup()
389  {  {
# Line 466  void DataLazy::LazyNodeSetup() Line 405  void DataLazy::LazyNodeSetup()
405    
406  // Creates an identity node  // Creates an identity node
407  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
408      : parent(p->getFunctionSpace(),p->getShape())          : parent(p->getFunctionSpace(),p->getShape())
409      ,m_sampleids(0),          ,m_sampleids(0),
410      m_samples(1)          m_samples(1)
411  {  {
412     if (p->isLazy())     if (p->isLazy())
413     {     {
414      // I don't want identity of Lazy.          // I don't want identity of Lazy.
415      // Question: Why would that be so bad?          // Question: Why would that be so bad?
416      // 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
417      throw DataException("Programmer error - attempt to create identity from a DataLazy.");          throw DataException("Programmer error - attempt to create identity from a DataLazy.");
418     }     }
419     else     else
420     {     {
421      p->makeLazyShared();          DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
422      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);          makeIdentity(dr);
     makeIdentity(dr);  
423  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
424     }     }
425  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
426  }  }
427    
428  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
429      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),          : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
430      m_op(op),          m_op(op),
431      m_axis_offset(0),          m_axis_offset(0),
432      m_transpose(0),          m_transpose(0),
433      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
434  {  {
435     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
436     {     {
437      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.");
438     }     }
439    
440     DataLazy_ptr lleft;     DataLazy_ptr lleft;
441     if (!left->isLazy())     if (!left->isLazy())
442     {     {
443      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
444     }     }
445     else     else
446     {     {
447      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
448     }     }
449     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
450     m_left=lleft;     m_left=lleft;
# Line 520  DataLazy::DataLazy(DataAbstract_ptr left Line 458  DataLazy::DataLazy(DataAbstract_ptr left
458    
459  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
460  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
461      : parent(resultFS(left,right,op), resultShape(left,right,op)),          : parent(resultFS(left,right,op), resultShape(left,right,op)),
462      m_op(op),          m_op(op),
463      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
464  {  {
465  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
466     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
467     {     {
468      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.");
469     }     }
470    
471     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
472     {     {
473      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
474      Data ltemp(left);          Data ltemp(left);
475      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
476      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
477     }     }
478     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
479     {     {
480      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
481      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
482  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
483     }     }
484     left->operandCheck(*right);     left->operandCheck(*right);
485    
486     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
487     {     {
488      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
489  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
490     }     }
491     else     else
492     {     {
493      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
494  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
495     }     }
496     if (right->isLazy())     if (right->isLazy())
497     {     {
498      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
499  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
500     }     }
501     else     else
502     {     {
503      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
504  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
505     }     }
506     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
507     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
508     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
509     {     {
510      m_readytype='E';          m_readytype='E';
511     }     }
512     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
513     {     {
514      m_readytype='T';          m_readytype='T';
515     }     }
516     else     else
517     {     {
518      m_readytype='C';          m_readytype='C';
519     }     }
520     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
521     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
# Line 588  LAZYDEBUG(cout << "(3)Lazy created with Line 526  LAZYDEBUG(cout << "(3)Lazy created with
526  }  }
527    
528  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)
529      : 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)),
530      m_op(op),          m_op(op),
531      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
532      m_transpose(transpose)          m_transpose(transpose)
533  {  {
534     if ((getOpgroup(op)!=G_TENSORPROD))     if ((getOpgroup(op)!=G_TENSORPROD))
535     {     {
536      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.");
537     }     }
538     if ((transpose>2) || (transpose<0))     if ((transpose>2) || (transpose<0))
539     {     {
540      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");
541     }     }
542     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
543     {     {
544      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
545      Data ltemp(left);          Data ltemp(left);
546      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
547      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
548     }     }
549     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
550     {     {
551      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
552      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
553     }     }
554  //    left->operandCheck(*right);  //    left->operandCheck(*right);
555    
556     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
557     {     {
558      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
559     }     }
560     else     else
561     {     {
562      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
563     }     }
564     if (right->isLazy())     if (right->isLazy())
565     {     {
566      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
567     }     }
568     else     else
569     {     {
570      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
571     }     }
572     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
573     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
574     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
575     {     {
576      m_readytype='E';          m_readytype='E';
577     }     }
578     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
579     {     {
580      m_readytype='T';          m_readytype='T';
581     }     }
582     else     else
583     {     {
584      m_readytype='C';          m_readytype='C';
585     }     }
586     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
587     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
# Line 655  LAZYDEBUG(cout << "(4)Lazy created with Line 593  LAZYDEBUG(cout << "(4)Lazy created with
593    
594    
595  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
596      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),          : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
597      m_op(op),          m_op(op),
598      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
599      m_transpose(0),          m_transpose(0),
600      m_tol(0)          m_tol(0)
601  {  {
602     if ((getOpgroup(op)!=G_NP1OUT_P))     if ((getOpgroup(op)!=G_NP1OUT_P))
603     {     {
604      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.");
605     }     }
606     DataLazy_ptr lleft;     DataLazy_ptr lleft;
607     if (!left->isLazy())     if (!left->isLazy())
608     {     {
609      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
610     }     }
611     else     else
612     {     {
613      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
614     }     }
615     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
616     m_left=lleft;     m_left=lleft;
# Line 685  LAZYDEBUG(cout << "(5)Lazy created with Line 623  LAZYDEBUG(cout << "(5)Lazy created with
623  }  }
624    
625  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
626      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
627      m_op(op),          m_op(op),
628      m_axis_offset(0),          m_axis_offset(0),
629      m_transpose(0),          m_transpose(0),
630      m_tol(tol)          m_tol(tol)
631  {  {
632     if ((getOpgroup(op)!=G_UNARY_P))     if ((getOpgroup(op)!=G_UNARY_P))
633     {     {
634      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.");
635     }     }
636     DataLazy_ptr lleft;     DataLazy_ptr lleft;
637     if (!left->isLazy())     if (!left->isLazy())
638     {     {
639      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
640     }     }
641     else     else
642     {     {
643      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
644     }     }
645     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
646     m_left=lleft;     m_left=lleft;
# Line 716  LAZYDEBUG(cout << "(6)Lazy created with Line 654  LAZYDEBUG(cout << "(6)Lazy created with
654    
655    
656  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)
657      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),          : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
658      m_op(op),          m_op(op),
659      m_axis_offset(axis0),          m_axis_offset(axis0),
660      m_transpose(axis1),          m_transpose(axis1),
661      m_tol(0)          m_tol(0)
662  {  {
663     if ((getOpgroup(op)!=G_NP1OUT_2P))     if ((getOpgroup(op)!=G_NP1OUT_2P))
664     {     {
665      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.");
666     }     }
667     DataLazy_ptr lleft;     DataLazy_ptr lleft;
668     if (!left->isLazy())     if (!left->isLazy())
669     {     {
670      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
671     }     }
672     else     else
673     {     {
674      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
675     }     }
676     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
677     m_left=lleft;     m_left=lleft;
# Line 751  namespace Line 689  namespace
689    
690      inline int max3(int a, int b, int c)      inline int max3(int a, int b, int c)
691      {      {
692      int t=(a>b?a:b);          int t=(a>b?a:b);
693      return (t>c?t:c);          return (t>c?t:c);
694    
695      }      }
696  }  }
697    
698  DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)  DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
699      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
700      m_op(CONDEVAL),          m_op(CONDEVAL),
701      m_axis_offset(0),          m_axis_offset(0),
702      m_transpose(0),          m_transpose(0),
703      m_tol(0)          m_tol(0)
704  {  {
705    
706     DataLazy_ptr lmask;     DataLazy_ptr lmask;
# Line 770  DataLazy::DataLazy(DataAbstract_ptr mask Line 708  DataLazy::DataLazy(DataAbstract_ptr mask
708     DataLazy_ptr lright;     DataLazy_ptr lright;
709     if (!mask->isLazy())     if (!mask->isLazy())
710     {     {
711      lmask=DataLazy_ptr(new DataLazy(mask));          lmask=DataLazy_ptr(new DataLazy(mask));
712     }     }
713     else     else
714     {     {
715      lmask=dynamic_pointer_cast<DataLazy>(mask);          lmask=dynamic_pointer_cast<DataLazy>(mask);
716     }     }
717     if (!left->isLazy())     if (!left->isLazy())
718     {     {
719      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
720     }     }
721     else     else
722     {     {
723      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
724     }     }
725     if (!right->isLazy())     if (!right->isLazy())
726     {     {
727      lright=DataLazy_ptr(new DataLazy(right));          lright=DataLazy_ptr(new DataLazy(right));
728     }     }
729     else     else
730     {     {
731      lright=dynamic_pointer_cast<DataLazy>(right);          lright=dynamic_pointer_cast<DataLazy>(right);
732     }     }
733     m_readytype=lmask->m_readytype;     m_readytype=lmask->m_readytype;
734     if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))     if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
735     {     {
736      throw DataException("Programmer Error - condEval arguments must have the same readytype");          throw DataException("Programmer Error - condEval arguments must have the same readytype");
737     }     }
738     m_left=lleft;     m_left=lleft;
739     m_right=lright;     m_right=lright;
# Line 822  DataLazy::~DataLazy() Line 760  DataLazy::~DataLazy()
760    For reasons of efficiency do not call this method on DataExpanded nodes.    For reasons of efficiency do not call this method on DataExpanded nodes.
761  */  */
762  DataReady_ptr  DataReady_ptr
763  DataLazy::collapseToReady()  DataLazy::collapseToReady() const
764  {  {
765    if (m_readytype=='E')    if (m_readytype=='E')
766    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
767      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
768    }    }
769    if (m_op==IDENTITY)    if (m_op==IDENTITY)
# Line 843  DataLazy::collapseToReady() Line 781  DataLazy::collapseToReady()
781    switch(m_op)    switch(m_op)
782    {    {
783      case ADD:      case ADD:
784      result=left+right;          result=left+right;
785      break;          break;
786      case SUB:            case SUB:          
787      result=left-right;          result=left-right;
788      break;          break;
789      case MUL:            case MUL:          
790      result=left*right;          result=left*right;
791      break;          break;
792      case DIV:            case DIV:          
793      result=left/right;          result=left/right;
794      break;          break;
795      case SIN:      case SIN:
796      result=left.sin();            result=left.sin();      
797      break;          break;
798      case COS:      case COS:
799      result=left.cos();          result=left.cos();
800      break;          break;
801      case TAN:      case TAN:
802      result=left.tan();          result=left.tan();
803      break;          break;
804      case ASIN:      case ASIN:
805      result=left.asin();          result=left.asin();
806      break;          break;
807      case ACOS:      case ACOS:
808      result=left.acos();          result=left.acos();
809      break;          break;
810      case ATAN:      case ATAN:
811      result=left.atan();          result=left.atan();
812      break;          break;
813      case SINH:      case SINH:
814      result=left.sinh();          result=left.sinh();
815      break;          break;
816      case COSH:      case COSH:
817      result=left.cosh();          result=left.cosh();
818      break;          break;
819      case TANH:      case TANH:
820      result=left.tanh();          result=left.tanh();
821      break;          break;
822      case ERF:      case ERF:
823      result=left.erf();          result=left.erf();
824      break;          break;
825     case ASINH:     case ASINH:
826      result=left.asinh();          result=left.asinh();
827      break;          break;
828     case ACOSH:     case ACOSH:
829      result=left.acosh();          result=left.acosh();
830      break;          break;
831     case ATANH:     case ATANH:
832      result=left.atanh();          result=left.atanh();
833      break;          break;
834      case LOG10:      case LOG10:
835      result=left.log10();          result=left.log10();
836      break;          break;
837      case LOG:      case LOG:
838      result=left.log();          result=left.log();
839      break;          break;
840      case SIGN:      case SIGN:
841      result=left.sign();          result=left.sign();
842      break;          break;
843      case ABS:      case ABS:
844      result=left.abs();          result=left.abs();
845      break;          break;
846      case NEG:      case NEG:
847      result=left.neg();          result=left.neg();
848      break;          break;
849      case POS:      case POS:
850      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
851      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
852      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
853      break;          break;
854      case EXP:      case EXP:
855      result=left.exp();          result=left.exp();
856      break;          break;
857      case SQRT:      case SQRT:
858      result=left.sqrt();          result=left.sqrt();
859      break;          break;
860      case RECIP:      case RECIP:
861      result=left.oneOver();          result=left.oneOver();
862      break;          break;
863      case GZ:      case GZ:
864      result=left.wherePositive();          result=left.wherePositive();
865      break;          break;
866      case LZ:      case LZ:
867      result=left.whereNegative();          result=left.whereNegative();
868      break;          break;
869      case GEZ:      case GEZ:
870      result=left.whereNonNegative();          result=left.whereNonNegative();
871      break;          break;
872      case LEZ:      case LEZ:
873      result=left.whereNonPositive();          result=left.whereNonPositive();
874      break;          break;
875      case NEZ:      case NEZ:
876      result=left.whereNonZero(m_tol);          result=left.whereNonZero(m_tol);
877      break;          break;
878      case EZ:      case EZ:
879      result=left.whereZero(m_tol);          result=left.whereZero(m_tol);
880      break;          break;
881      case SYM:      case SYM:
882      result=left.symmetric();          result=left.symmetric();
883      break;          break;
884      case NSYM:      case NSYM:
885      result=left.nonsymmetric();          result=left.antisymmetric();
886      break;          break;
887      case PROD:      case PROD:
888      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);          result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
889      break;          break;
890      case TRANS:      case TRANS:
891      result=left.transpose(m_axis_offset);          result=left.transpose(m_axis_offset);
892      break;          break;
893      case TRACE:      case TRACE:
894      result=left.trace(m_axis_offset);          result=left.trace(m_axis_offset);
895      break;          break;
896      case SWAP:      case SWAP:
897      result=left.swapaxes(m_axis_offset, m_transpose);          result=left.swapaxes(m_axis_offset, m_transpose);
898      break;          break;
899      case MINVAL:      case MINVAL:
900      result=left.minval();          result=left.minval();
901      break;          break;
902      case MAXVAL:      case MAXVAL:
903      result=left.minval();          result=left.minval();
904            break;
905        case HER:
906        result=left.hermitian();
907      break;      break;
908      default:      default:
909      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)+".");
910    }    }
911    return result.borrowReadyPtr();    return result.borrowReadyPtr();
912  }  }
# Line 977  DataLazy::collapseToReady() Line 918  DataLazy::collapseToReady()
918     the purpose of using DataLazy in the first place).     the purpose of using DataLazy in the first place).
919  */  */
920  void  void
921  DataLazy::collapse()  DataLazy::collapse() const
922  {  {
923    if (m_op==IDENTITY)    if (m_op==IDENTITY)
924    {    {
925      return;          return;
926    }    }
927    if (m_readytype=='E')    if (m_readytype=='E')
928    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
929      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
930    }    }
931    m_id=collapseToReady();    m_id=collapseToReady();
932    m_op=IDENTITY;    m_op=IDENTITY;
933  }  }
934    
   
   
   
   
   
 #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;\  
     }  
   
   
935  // The result will be stored in m_samples  // The result will be stored in m_samples
936  // 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
937  const DataTypes::ValueType*  const DataTypes::RealVectorType*
938  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
939  {  {
940  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
941      // 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
942    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
943    {    {
944      collapse();          collapse();
945    }    }
946    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
947    {    {
948      const ValueType& vec=m_id->getVectorRO();      const RealVectorType& vec=m_id->getVectorRO();
949      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
950  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
951  int x;  int x;
# Line 1043  if (&x<stackend[omp_get_thread_num()]) Line 962  if (&x<stackend[omp_get_thread_num()])
962    }    }
963    if (m_sampleids[tid]==sampleNo)    if (m_sampleids[tid]==sampleNo)
964    {    {
965      roffset=tid*m_samplesize;          roffset=tid*m_samplesize;
966      return &(m_samples);        // sample is already resolved          return &(m_samples);            // sample is already resolved
967    }    }
968    m_sampleids[tid]=sampleNo;    m_sampleids[tid]=sampleNo;
969    
# Line 1064  if (&x<stackend[omp_get_thread_num()]) Line 983  if (&x<stackend[omp_get_thread_num()])
983    }    }
984  }  }
985    
986  const DataTypes::ValueType*  const DataTypes::RealVectorType*
987  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
988  {  {
989      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
990      // 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
991      // processing single points.          // processing single points.
992      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
993    if (m_readytype!='E')    if (m_readytype!='E')
994    {    {
995      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 1079  DataLazy::resolveNodeUnary(int tid, int Line 998  DataLazy::resolveNodeUnary(int tid, int
998    {    {
999      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1000    }    }
1001    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1002    const double* left=&((*leftres)[roffset]);    const double* left=&((*leftres)[roffset]);
1003    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1004    double* result=&(m_samples[roffset]);    double* result=&(m_samples[roffset]);
1005    switch (m_op)    if (m_op==POS)
1006    {    {
1007      case SIN:        // this should be prevented earlier
1008      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);      // operation is meaningless for lazy
1009      break;          throw DataException("Programmer error - POS not supported for lazy data.");    
1010      case COS:    }
1011      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);    tensor_unary_array_operation(m_samplesize,
1012      break;                               left,
1013      case TAN:                               result,
1014      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);                               m_op,
1015      break;                               m_tol);  
     case ASIN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);  
     break;  
     case ACOS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);  
     break;  
     case ATAN:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);  
     break;  
     case SINH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);  
     break;  
     case COSH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);  
     break;  
     case TANH:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);  
     break;  
     case ERF:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::erf);  
     break;  
 #endif  
    case ASINH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 #endif    
     break;  
    case ACOSH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 #endif    
     break;  
    case ATANH:  
 #if defined (_WIN32) && !defined(__INTEL_COMPILER)  
     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 #endif    
     break;  
     case LOG10:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);  
     break;  
     case LOG:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);  
     break;  
     case SIGN:  
     tensor_unary_operation(m_samplesize, left, result, escript::fsign);  
     break;  
     case ABS:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);  
     break;  
     case NEG:  
     tensor_unary_operation(m_samplesize, left, result, negate<double>());  
     break;  
     case POS:  
     // it doesn't mean anything for delayed.  
     // it will just trigger a deep copy of the lazy object  
     throw DataException("Programmer error - POS not supported for lazy data.");  
     break;  
     case EXP:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);  
     break;  
     case SQRT:  
     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);  
     break;  
     case RECIP:  
     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));  
     break;  
     case GZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));  
     break;  
     case LZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));  
     break;  
     case GEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));  
     break;  
     case LEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));  
     break;  
 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently  
     case NEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));  
     break;  
     case EZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));  
     break;  
   
     default:  
     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");  
   }  
1016    return &(m_samples);    return &(m_samples);
1017  }  }
1018    
1019    
1020  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1021  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1022  {  {
1023      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1024      // 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
1025      // processing single points.          // processing single points.
1026      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1027    if (m_readytype!='E')    if (m_readytype!='E')
1028    {    {
1029      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 1212  DataLazy::resolveNodeReduction(int tid, Line 1033  DataLazy::resolveNodeReduction(int tid,
1033      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1034    }    }
1035    size_t loffset=0;    size_t loffset=0;
1036    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1037    
1038    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1039    unsigned int ndpps=getNumDPPSample();    unsigned int ndpps=getNumDPPSample();
# Line 1221  DataLazy::resolveNodeReduction(int tid, Line 1042  DataLazy::resolveNodeReduction(int tid,
1042    switch (m_op)    switch (m_op)
1043    {    {
1044      case MINVAL:      case MINVAL:
1045      {          {
1046        for (unsigned int z=0;z<ndpps;++z)            for (unsigned int z=0;z<ndpps;++z)
1047        {            {
1048          FMin op;              FMin op;
1049          *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());              *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1050          loffset+=psize;              loffset+=psize;
1051          result++;              result++;
1052        }            }
1053      }          }
1054      break;          break;
1055      case MAXVAL:      case MAXVAL:
1056      {          {
1057        for (unsigned int z=0;z<ndpps;++z)            for (unsigned int z=0;z<ndpps;++z)
1058        {            {
1059        FMax op;            FMax op;
1060        *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);            *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1061        loffset+=psize;            loffset+=psize;
1062        result++;            result++;
1063        }            }
1064      }          }
1065      break;          break;
1066      default:      default:
1067      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1068    }    }
1069    return &(m_samples);    return &(m_samples);
1070  }  }
1071    
1072  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1073  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1074  {  {
1075      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1076      // 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
1077      // processing single points.          // processing single points.
1078    if (m_readytype!='E')    if (m_readytype!='E')
1079    {    {
1080      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 1263  DataLazy::resolveNodeNP1OUT(int tid, int Line 1084  DataLazy::resolveNodeNP1OUT(int tid, int
1084      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1085    }    }
1086    size_t subroffset;    size_t subroffset;
1087    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1088    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1089    size_t loop=0;    size_t loop=0;
1090    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
# Line 1272  DataLazy::resolveNodeNP1OUT(int tid, int Line 1093  DataLazy::resolveNodeNP1OUT(int tid, int
1093    switch (m_op)    switch (m_op)
1094    {    {
1095      case SYM:      case SYM:
1096      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1097      {          {
1098          DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              escript::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1099          subroffset+=step;              subroffset+=step;
1100          offset+=step;              offset+=step;
1101      }          }
1102      break;          break;
1103      case NSYM:      case NSYM:
1104      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1105      {          {
1106          DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              escript::antisymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1107          subroffset+=step;              subroffset+=step;
1108          offset+=step;              offset+=step;
1109      }          }
1110      break;          break;
1111      default:      default:
1112      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1113    }    }
1114    return &m_samples;    return &m_samples;
1115  }  }
1116    
1117  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1118  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1119  {  {
1120      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1121      // 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
1122      // processing single points.          // processing single points.
1123    if (m_readytype!='E')    if (m_readytype!='E')
1124    {    {
1125      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 1309  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1130  DataLazy::resolveNodeNP1OUT_P(int tid, i
1130    }    }
1131    size_t subroffset;    size_t subroffset;
1132    size_t offset;    size_t offset;
1133    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1134    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1135    offset=roffset;    offset=roffset;
1136    size_t loop=0;    size_t loop=0;
# Line 1319  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1140  DataLazy::resolveNodeNP1OUT_P(int tid, i
1140    switch (m_op)    switch (m_op)
1141    {    {
1142      case TRACE:      case TRACE:
1143      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1144      {          {
1145              DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);              escript::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1146          subroffset+=instep;              subroffset+=instep;
1147          offset+=outstep;              offset+=outstep;
1148      }          }
1149      break;          break;
1150      case TRANS:      case TRANS:
1151      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1152      {          {
1153              DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);              escript::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1154          subroffset+=instep;              subroffset+=instep;
1155          offset+=outstep;              offset+=outstep;
1156      }          }
1157      break;          break;
1158      default:      default:
1159      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1160    }    }
1161    return &m_samples;    return &m_samples;
1162  }  }
1163    
1164    
1165  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1166  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1167  {  {
1168    if (m_readytype!='E')    if (m_readytype!='E')
1169    {    {
# Line 1354  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1175  DataLazy::resolveNodeNP1OUT_2P(int tid,
1175    }    }
1176    size_t subroffset;    size_t subroffset;
1177    size_t offset;    size_t offset;
1178    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1179    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1180    offset=roffset;    offset=roffset;
1181    size_t loop=0;    size_t loop=0;
# Line 1364  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1185  DataLazy::resolveNodeNP1OUT_2P(int tid,
1185    switch (m_op)    switch (m_op)
1186    {    {
1187      case SWAP:      case SWAP:
1188      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1189      {          {
1190              DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);              escript::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1191          subroffset+=instep;              subroffset+=instep;
1192          offset+=outstep;              offset+=outstep;
1193      }          }
1194      break;          break;
1195      default:      default:
1196      throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1197    }    }
1198    return &m_samples;    return &m_samples;
1199  }  }
1200    
1201  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1202  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1203  {  {
1204    if (m_readytype!='E')    if (m_readytype!='E')
1205    {    {
# Line 1390  DataLazy::resolveNodeCondEval(int tid, i Line 1211  DataLazy::resolveNodeCondEval(int tid, i
1211    }    }
1212    size_t subroffset;    size_t subroffset;
1213    
1214    const ValueType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1215    const ValueType* srcres=0;    const RealVectorType* srcres=0;
1216    if ((*maskres)[subroffset]>0)    if ((*maskres)[subroffset]>0)
1217    {    {
1218      srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1219    }    }
1220    else    else
1221    {    {
1222      srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1223    }    }
1224    
1225    // Now we need to copy the result    // Now we need to copy the result
# Line 1406  DataLazy::resolveNodeCondEval(int tid, i Line 1227  DataLazy::resolveNodeCondEval(int tid, i
1227    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1228    for (int i=0;i<m_samplesize;++i)    for (int i=0;i<m_samplesize;++i)
1229    {    {
1230      m_samples[roffset+i]=(*srcres)[subroffset+i];            m_samples[roffset+i]=(*srcres)[subroffset+i];  
1231    }    }
1232    
1233    return &m_samples;    return &m_samples;
# Line 1421  DataLazy::resolveNodeCondEval(int tid, i Line 1242  DataLazy::resolveNodeCondEval(int tid, i
1242  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1243  // For example, 2+Vector.  // For example, 2+Vector.
1244  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1245  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1246  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1247  {  {
1248  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1249    
1250    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
1251      // first work out which of the children are expanded          // first work out which of the children are expanded
1252    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1253    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1254    if (!leftExp && !rightExp)    if (!leftExp && !rightExp)
1255    {    {
1256      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'.");
1257    }    }
1258    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1259    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1260    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1261    {    {
1262      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.");
1263    }    }
1264    size_t leftsize=m_left->getNoValues();    size_t leftsize=m_left->getNoValues();
1265    size_t rightsize=m_right->getNoValues();    size_t rightsize=m_right->getNoValues();
1266    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
1267    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
1268    int rightstep=0;    int rightstep=0;
1269    int numsteps=0;       // total number of steps for the inner loop    int numsteps=0;               // total number of steps for the inner loop
1270    int oleftstep=0;  // the o variables refer to the outer loop    int oleftstep=0;      // the o variables refer to the outer loop
1271    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
1272    int onumsteps=1;    int onumsteps=1;
1273        
1274    bool LES=(leftExp && leftScalar); // Left is an expanded scalar    bool LES=(leftExp && leftScalar);     // Left is an expanded scalar
1275    bool RES=(rightExp && rightScalar);    bool RES=(rightExp && rightScalar);
1276    bool LS=(!leftExp && leftScalar); // left is a single scalar    bool LS=(!leftExp && leftScalar);     // left is a single scalar
1277    bool RS=(!rightExp && rightScalar);    bool RS=(!rightExp && rightScalar);
1278    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1279    bool RN=(!rightExp && !rightScalar);    bool RN=(!rightExp && !rightScalar);
1280    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1281    bool REN=(rightExp && !rightScalar);    bool REN=(rightExp && !rightScalar);
1282    
1283    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
1284    {    {
1285      chunksize=m_left->getNumDPPSample()*leftsize;          chunksize=m_left->getNumDPPSample()*leftsize;
1286      leftstep=0;          leftstep=0;
1287      rightstep=0;          rightstep=0;
1288      numsteps=1;          numsteps=1;
1289    }    }
1290    else if (LES || RES)    else if (LES || RES)
1291    {    {
1292      chunksize=1;          chunksize=1;
1293      if (LES)        // left is an expanded scalar          if (LES)                // left is an expanded scalar
1294      {          {
1295          if (RS)                  if (RS)
1296          {                  {
1297             leftstep=1;                     leftstep=1;
1298             rightstep=0;                     rightstep=0;
1299             numsteps=m_left->getNumDPPSample();                     numsteps=m_left->getNumDPPSample();
1300          }                  }
1301          else        // RN or REN                  else            // RN or REN
1302          {                  {
1303             leftstep=0;                     leftstep=0;
1304             oleftstep=1;                     oleftstep=1;
1305             rightstep=1;                     rightstep=1;
1306             orightstep=(RN ? -(int)rightsize : 0);                     orightstep=(RN ? -(int)rightsize : 0);
1307             numsteps=rightsize;                     numsteps=rightsize;
1308             onumsteps=m_left->getNumDPPSample();                     onumsteps=m_left->getNumDPPSample();
1309          }                  }
1310      }          }
1311      else        // right is an expanded scalar          else            // right is an expanded scalar
1312      {          {
1313          if (LS)                  if (LS)
1314          {                  {
1315             rightstep=1;                     rightstep=1;
1316             leftstep=0;                     leftstep=0;
1317             numsteps=m_right->getNumDPPSample();                     numsteps=m_right->getNumDPPSample();
1318          }                  }
1319          else                  else
1320          {                  {
1321             rightstep=0;                     rightstep=0;
1322             orightstep=1;                     orightstep=1;
1323             leftstep=1;                     leftstep=1;
1324             oleftstep=(LN ? -(int)leftsize : 0);                     oleftstep=(LN ? -(int)leftsize : 0);
1325             numsteps=leftsize;                     numsteps=leftsize;
1326             onumsteps=m_right->getNumDPPSample();                     onumsteps=m_right->getNumDPPSample();
1327          }                  }
1328      }          }
1329    }    }
1330    else  // this leaves (LEN, RS), (LEN, RN) and their transposes    else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1331    {    {
1332      if (LEN)    // and Right will be a single value          if (LEN)        // and Right will be a single value
1333      {          {
1334          chunksize=rightsize;                  chunksize=rightsize;
1335          leftstep=rightsize;                  leftstep=rightsize;
1336          rightstep=0;                  rightstep=0;
1337          numsteps=m_left->getNumDPPSample();                  numsteps=m_left->getNumDPPSample();
1338          if (RS)                  if (RS)
1339          {                  {
1340             numsteps*=leftsize;                     numsteps*=leftsize;
1341          }                  }
1342      }          }
1343      else    // REN          else    // REN
1344      {          {
1345          chunksize=leftsize;                  chunksize=leftsize;
1346          rightstep=leftsize;                  rightstep=leftsize;
1347          leftstep=0;                  leftstep=0;
1348          numsteps=m_right->getNumDPPSample();                  numsteps=m_right->getNumDPPSample();
1349          if (LS)                  if (LS)
1350          {                  {
1351             numsteps*=rightsize;                     numsteps*=rightsize;
1352          }                  }
1353      }          }
1354    }    }
1355    
1356    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1357      // Get the values of sub-expressions          // Get the values of sub-expressions
1358    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1359    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1360  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1361  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;)
1362  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
# Line 1549  LAZYDEBUG(cout << "Right res["<< rroffse Line 1370  LAZYDEBUG(cout << "Right res["<< rroffse
1370    
1371    
1372    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1373    double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we received    double* resultp=&(m_samples[roffset]);                // results are stored at the vector offset we received
1374    switch(m_op)    switch(m_op)
1375    {    {
1376      case ADD:      case ADD:
1377          PROC_OP(NO_ARG,plus<double>());          //PROC_OP(NO_ARG,plus<double>());
1378      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1379                 &(*left)[0],
1380                 &(*right)[0],
1381                 chunksize,
1382                 onumsteps,
1383                 numsteps,
1384                 resultStep,
1385                 leftstep,
1386                 rightstep,
1387                 oleftstep,
1388                 orightstep,
1389                 lroffset,
1390                 rroffset,
1391                 escript::ES_optype::ADD);  
1392            break;
1393      case SUB:      case SUB:
1394      PROC_OP(NO_ARG,minus<double>());        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1395      break;               &(*left)[0],
1396                 &(*right)[0],
1397                 chunksize,
1398                 onumsteps,
1399                 numsteps,
1400                 resultStep,
1401                 leftstep,
1402                 rightstep,
1403                 oleftstep,
1404                 orightstep,
1405                 lroffset,
1406                 rroffset,
1407                 escript::ES_optype::SUB);        
1408            //PROC_OP(NO_ARG,minus<double>());
1409            break;
1410      case MUL:      case MUL:
1411      PROC_OP(NO_ARG,multiplies<double>());          //PROC_OP(NO_ARG,multiplies<double>());
1412      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1413                 &(*left)[0],
1414                 &(*right)[0],
1415                 chunksize,
1416                 onumsteps,
1417                 numsteps,
1418                 resultStep,
1419                 leftstep,
1420                 rightstep,
1421                 oleftstep,
1422                 orightstep,
1423                 lroffset,
1424                 rroffset,
1425                 escript::ES_optype::MUL);        
1426            break;
1427      case DIV:      case DIV:
1428      PROC_OP(NO_ARG,divides<double>());          //PROC_OP(NO_ARG,divides<double>());
1429      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1430                 &(*left)[0],
1431                 &(*right)[0],
1432                 chunksize,
1433                 onumsteps,
1434                 numsteps,
1435                 resultStep,
1436                 leftstep,
1437                 rightstep,
1438                 oleftstep,
1439                 orightstep,
1440                 lroffset,
1441                 rroffset,
1442                 escript::ES_optype::DIV);        
1443            break;
1444      case POW:      case POW:
1445         PROC_OP(double (double,double),::pow);         //PROC_OP(double (double,double),::pow);
1446      break;        escript::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1447                 &(*left)[0],
1448                 &(*right)[0],
1449                 chunksize,
1450                 onumsteps,
1451                 numsteps,
1452                 resultStep,
1453                 leftstep,
1454                 rightstep,
1455                 oleftstep,
1456                 orightstep,
1457                 lroffset,
1458                 rroffset,
1459                 escript::ES_optype::POW);        
1460            break;
1461      default:      default:
1462      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1463    }    }
1464  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1465    return &m_samples;    return &m_samples;
# Line 1578  LAZYDEBUG(cout << "Result res[" << roffs Line 1469  LAZYDEBUG(cout << "Result res[" << roffs
1469  // 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
1470  // 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.
1471  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1472  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1473  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1474  {  {
1475  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1476    
1477    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
1478      // first work out which of the children are expanded          // first work out which of the children are expanded
1479    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1480    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1481    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1482    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
1483    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1484    
1485    int resultStep=getNoValues();    int resultStep=getNoValues();
1486    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1487    size_t offset=roffset;    size_t offset=roffset;
1488    
1489    const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);    const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1490    
1491    const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1492    
1493  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;
1494  cout << getNoValues() << endl;)  cout << getNoValues() << endl;)
# Line 1611  LAZYDEBUG(cout << "m_samplesize=" << m_s Line 1502  LAZYDEBUG(cout << "m_samplesize=" << m_s
1502  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1503  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1504    
1505    double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we received    double* resultp=&(m_samples[offset]);         // results are stored at the vector offset we received
1506    switch(m_op)    switch(m_op)
1507    {    {
1508      case PROD:      case PROD:
1509      for (int i=0;i<steps;++i,resultp+=resultStep)          for (int i=0;i<steps;++i,resultp+=resultStep)
1510      {          {
1511            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1512            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1513    
1514  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1515  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1516    
1517            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);
1518    
1519        lroffset+=leftStep;            lroffset+=leftStep;
1520        rroffset+=rightStep;            rroffset+=rightStep;
1521      }          }
1522      break;          break;
1523      default:      default:
1524      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1525    }    }
1526    roffset=offset;    roffset=offset;
1527    return &m_samples;    return &m_samples;
1528  }  }
1529    
1530    
1531  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1532  DataLazy::resolveSample(int sampleNo, size_t& roffset)  DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1533  {  {
1534  #ifdef _OPENMP  #ifdef _OPENMP
1535      int tid=omp_get_thread_num();          int tid=omp_get_thread_num();
1536  #else  #else
1537      int tid=0;          int tid=0;
1538  #endif  #endif
1539    
1540  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1541      stackstart[tid]=&tid;          stackstart[tid]=&tid;
1542      stackend[tid]=&tid;          stackend[tid]=&tid;
1543      const DataTypes::ValueType* r=resolveNodeSample(tid, sampleNo, roffset);          const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1544      size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];          size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1545      #pragma omp critical          #pragma omp critical
1546      if (d>maxstackuse)          if (d>maxstackuse)
1547      {          {
1548  cout << "Max resolve Stack use " << d << endl;  cout << "Max resolve Stack use " << d << endl;
1549          maxstackuse=d;                  maxstackuse=d;
1550      }          }
1551      return r;          return r;
1552  #else  #else
1553      return resolveNodeSample(tid, sampleNo, roffset);          return resolveNodeSample(tid, sampleNo, roffset);
1554  #endif  #endif
1555  }  }
1556    
# Line 1669  void Line 1560  void
1560  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
1561  {  {
1562     if (m_op==IDENTITY)     if (m_op==IDENTITY)
1563      return;          return;
1564     DataReady_ptr p=resolveNodeWorker();     DataReady_ptr p=resolveNodeWorker();
1565     makeIdentity(p);     makeIdentity(p);
1566  }  }
# Line 1706  DataLazy::resolveGroupWorker(std::vector Line 1597  DataLazy::resolveGroupWorker(std::vector
1597  {  {
1598    if (dats.empty())    if (dats.empty())
1599    {    {
1600      return;          return;
1601    }    }
1602    vector<DataLazy*> work;    vector<DataLazy*> work;
1603    FunctionSpace fs=dats[0]->getFunctionSpace();    FunctionSpace fs=dats[0]->getFunctionSpace();
1604    bool match=true;    bool match=true;
1605    for (int i=dats.size()-1;i>=0;--i)    for (int i=dats.size()-1;i>=0;--i)
1606    {    {
1607      if (dats[i]->m_readytype!='E')          if (dats[i]->m_readytype!='E')
1608      {          {
1609          dats[i]->collapse();                  dats[i]->collapse();
1610      }          }
1611      if (dats[i]->m_op!=IDENTITY)          if (dats[i]->m_op!=IDENTITY)
1612      {          {
1613          work.push_back(dats[i]);                  work.push_back(dats[i]);
1614          if (fs!=dats[i]->getFunctionSpace())                  if (fs!=dats[i]->getFunctionSpace())
1615          {                  {
1616              match=false;                          match=false;
1617          }                  }
1618      }          }
1619    }    }
1620    if (work.empty())    if (work.empty())
1621    {    {
1622      return;     // no work to do          return;         // no work to do
1623    }    }
1624    if (match)    // all functionspaces match.  Yes I realise this is overly strict    if (match)    // all functionspaces match.  Yes I realise this is overly strict
1625    {     // it is possible that dats[0] is one of the objects which we discarded and    {             // it is possible that dats[0] is one of the objects which we discarded and
1626          // all the other functionspaces match.                  // all the other functionspaces match.
1627      vector<DataExpanded*> dep;          vector<DataExpanded*> dep;
1628      vector<ValueType*> vecs;          vector<RealVectorType*> vecs;
1629      for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1630      {          {
1631          dep.push_back(new DataExpanded(fs,work[i]->getShape(), ValueType(work[i]->getNoValues())));                  dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
1632          vecs.push_back(&(dep[i]->getVectorRW()));                  vecs.push_back(&(dep[i]->getVectorRW()));
1633      }          }
1634      int totalsamples=work[0]->getNumSamples();          int totalsamples=work[0]->getNumSamples();
1635      const ValueType* res=0; // Storage for answer          const RealVectorType* res=0; // Storage for answer
1636      int sample;          int sample;
1637      #pragma omp parallel private(sample, res)          #pragma omp parallel private(sample, res)
1638      {          {
1639          size_t roffset=0;              size_t roffset=0;
1640          #pragma omp for schedule(static)              #pragma omp for schedule(static)
1641          for (sample=0;sample<totalsamples;++sample)              for (sample=0;sample<totalsamples;++sample)
1642          {              {
1643          roffset=0;                  roffset=0;
1644          int j;                  int j;
1645          for (j=work.size()-1;j>=0;--j)                  for (j=work.size()-1;j>=0;--j)
1646          {                  {
1647  #ifdef _OPENMP  #ifdef _OPENMP
1648                  res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);                      res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1649  #else  #else
1650                  res=work[j]->resolveNodeSample(0,sample,roffset);                      res=work[j]->resolveNodeSample(0,sample,roffset);
1651  #endif  #endif
1652                  DataVector::size_type outoffset=dep[j]->getPointOffset(sample,0);                      RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1653                  memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(DataVector::ElementType));                      memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
1654          }                  }
1655          }              }
1656      }          }
1657      // Now we need to load the new results as identity ops into the lazy nodes          // Now we need to load the new results as identity ops into the lazy nodes
1658      for (int i=work.size()-1;i>=0;--i)          for (int i=work.size()-1;i>=0;--i)
1659      {          {
1660          work[i]->makeIdentity(boost::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));              work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1661      }          }
1662    }    }
1663    else  // functionspaces do not match    else  // functionspaces do not match
1664    {    {
1665      for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1666      {          {
1667          work[i]->resolveToIdentity();                  work[i]->resolveToIdentity();
1668      }          }
1669    }    }
1670  }  }
1671    
# Line 1784  DataLazy::resolveGroupWorker(std::vector Line 1675  DataLazy::resolveGroupWorker(std::vector
1675  DataReady_ptr  DataReady_ptr
1676  DataLazy::resolveNodeWorker()  DataLazy::resolveNodeWorker()
1677  {  {
1678    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally    if (m_readytype!='E')         // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1679    {    {
1680      collapse();      collapse();
1681    }    }
1682    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.
1683    {    {
1684      return m_id;      return m_id;
1685    }    }
1686      // 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'
1687    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  RealVectorType(getNoValues()));
1688    ValueType& resvec=result->getVectorRW();    RealVectorType& resvec=result->getVectorRW();
1689    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1690    
1691    int sample;    int sample;
1692    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1693    const ValueType* res=0;   // Storage for answer    const RealVectorType* res=0;       // Storage for answer
1694  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1695    #pragma omp parallel private(sample,res)    #pragma omp parallel private(sample,res)
1696    {    {
1697      size_t roffset=0;          size_t roffset=0;
1698  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1699      stackstart[omp_get_thread_num()]=&roffset;          stackstart[omp_get_thread_num()]=&roffset;
1700      stackend[omp_get_thread_num()]=&roffset;          stackend[omp_get_thread_num()]=&roffset;
1701  #endif  #endif
1702      #pragma omp for schedule(static)          #pragma omp for schedule(static)
1703      for (sample=0;sample<totalsamples;++sample)          for (sample=0;sample<totalsamples;++sample)
1704      {          {
1705          roffset=0;                  roffset=0;
1706  #ifdef _OPENMP  #ifdef _OPENMP
1707              res=resolveNodeSample(omp_get_thread_num(),sample,roffset);                  res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1708  #else  #else
1709              res=resolveNodeSample(0,sample,roffset);                  res=resolveNodeSample(0,sample,roffset);
1710  #endif  #endif
1711  LAZYDEBUG(cout << "Sample #" << sample << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1712  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1713              DataVector::size_type outoffset=result->getPointOffset(sample,0);                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1714              memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));                  memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1715      }          }
1716    }    }
1717  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1718    for (int i=0;i<getNumberOfThreads();++i)    for (int i=0;i<getNumberOfThreads();++i)
1719    {    {
1720      size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);          size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1721  //  cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;  //      cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1722      if (r>maxstackuse)          if (r>maxstackuse)
1723      {          {
1724          maxstackuse=r;                  maxstackuse=r;
1725      }          }
1726    }    }
1727    cout << "Max resolve Stack use=" << maxstackuse << endl;    cout << "Max resolve Stack use=" << maxstackuse << endl;
1728  #endif  #endif
# Line 1843  DataLazy::toString() const Line 1734  DataLazy::toString() const
1734  {  {
1735    ostringstream oss;    ostringstream oss;
1736    oss << "Lazy Data: [depth=" << m_height<< "] ";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1737    switch (escriptParams.getLAZY_STR_FMT())    switch (escriptParams.getLazyStrFmt())
1738    {    {
1739    case 1:   // tree format    case 1:       // tree format
1740      oss << endl;          oss << endl;
1741      intoTreeString(oss,"");          intoTreeString(oss,"");
1742      break;          break;
1743    case 2:   // just the depth    case 2:       // just the depth
1744      break;          break;
1745    default:    default:
1746      intoString(oss);          intoString(oss);
1747      break;          break;
1748    }    }
1749    return oss.str();    return oss.str();
1750  }  }
# Line 1866  DataLazy::intoString(ostringstream& oss) Line 1757  DataLazy::intoString(ostringstream& oss)
1757    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1758    {    {
1759    case G_IDENTITY:    case G_IDENTITY:
1760      if (m_id->isExpanded())          if (m_id->isExpanded())
1761      {          {
1762         oss << "E";             oss << "E";
1763      }          }
1764      else if (m_id->isTagged())          else if (m_id->isTagged())
1765      {          {
1766        oss << "T";            oss << "T";
1767      }          }
1768      else if (m_id->isConstant())          else if (m_id->isConstant())
1769      {          {
1770        oss << "C";            oss << "C";
1771      }          }
1772      else          else
1773      {          {
1774        oss << "?";            oss << "?";
1775      }          }
1776      oss << '@' << m_id.get();          oss << '@' << m_id.get();
1777      break;          break;
1778    case G_BINARY:    case G_BINARY:
1779      oss << '(';          oss << '(';
1780      m_left->intoString(oss);          m_left->intoString(oss);
1781      oss << ' ' << opToString(m_op) << ' ';          oss << ' ' << opToString(m_op) << ' ';
1782      m_right->intoString(oss);          m_right->intoString(oss);
1783      oss << ')';          oss << ')';
1784      break;          break;
1785    case G_UNARY:    case G_UNARY:
1786    case G_UNARY_P:    case G_UNARY_P:
1787    case G_NP1OUT:    case G_NP1OUT:
1788    case G_NP1OUT_P:    case G_NP1OUT_P:
1789    case G_REDUCTION:    case G_REDUCTION:
1790      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1791      m_left->intoString(oss);          m_left->intoString(oss);
1792      oss << ')';          oss << ')';
1793      break;          break;
1794    case G_TENSORPROD:    case G_TENSORPROD:
1795      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1796      m_left->intoString(oss);          m_left->intoString(oss);
1797      oss << ", ";          oss << ", ";
1798      m_right->intoString(oss);          m_right->intoString(oss);
1799      oss << ')';          oss << ')';
1800      break;          break;
1801    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1802      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1803      m_left->intoString(oss);          m_left->intoString(oss);
1804      oss << ", " << m_axis_offset << ", " << m_transpose;          oss << ", " << m_axis_offset << ", " << m_transpose;
1805      oss << ')';          oss << ')';
1806      break;          break;
1807    case G_CONDEVAL:    case G_CONDEVAL:
1808      oss << opToString(m_op)<< '(' ;          oss << opToString(m_op)<< '(' ;
1809      m_mask->intoString(oss);          m_mask->intoString(oss);
1810      oss << " ? ";          oss << " ? ";
1811      m_left->intoString(oss);          m_left->intoString(oss);
1812      oss << " : ";          oss << " : ";
1813      m_right->intoString(oss);          m_right->intoString(oss);
1814      oss << ')';          oss << ')';
1815      break;          break;
1816    default:    default:
1817      oss << "UNKNOWN";          oss << "UNKNOWN";
1818    }    }
1819  }  }
1820    
# Line 1935  DataLazy::intoTreeString(ostringstream& Line 1826  DataLazy::intoTreeString(ostringstream&
1826    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1827    {    {
1828    case G_IDENTITY:    case G_IDENTITY:
1829      if (m_id->isExpanded())          if (m_id->isExpanded())
1830      {          {
1831         oss << "E";             oss << "E";
1832      }          }
1833      else if (m_id->isTagged())          else if (m_id->isTagged())
1834      {          {
1835        oss << "T";            oss << "T";
1836      }          }
1837      else if (m_id->isConstant())          else if (m_id->isConstant())
1838      {          {
1839        oss << "C";            oss << "C";
1840      }          }
1841      else          else
1842      {          {
1843        oss << "?";            oss << "?";
1844      }          }
1845      oss << '@' << m_id.get() << endl;          oss << '@' << m_id.get() << endl;
1846      break;          break;
1847    case G_BINARY:    case G_BINARY:
1848      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1849      indent+='.';          indent+='.';
1850      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1851      m_right->intoTreeString(oss, indent);          m_right->intoTreeString(oss, indent);
1852      break;          break;
1853    case G_UNARY:    case G_UNARY:
1854    case G_UNARY_P:    case G_UNARY_P:
1855    case G_NP1OUT:    case G_NP1OUT:
1856    case G_NP1OUT_P:    case G_NP1OUT_P:
1857    case G_REDUCTION:    case G_REDUCTION:
1858      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1859      indent+='.';          indent+='.';
1860      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1861      break;          break;
1862    case G_TENSORPROD:    case G_TENSORPROD:
1863      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1864      indent+='.';          indent+='.';
1865      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1866      m_right->intoTreeString(oss, indent);          m_right->intoTreeString(oss, indent);
1867      break;          break;
1868    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1869      oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;          oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1870      indent+='.';          indent+='.';
1871      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1872      break;          break;
1873    default:    default:
1874      oss << "UNKNOWN";          oss << "UNKNOWN";
1875    }    }
1876  }  }
1877    
1878    
1879  DataAbstract*  DataAbstract*
1880  DataLazy::deepCopy()  DataLazy::deepCopy() const
1881  {  {
1882    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1883    {    {
1884    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1885    case G_UNARY:    case G_UNARY:
1886    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1887    case G_UNARY_P:   return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);    case G_UNARY_P:       return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1888    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);    case G_BINARY:        return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1889    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);
1890    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);
1891    case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);    case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
1892    case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);    case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1893    default:    default:
1894      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)+".");
1895    }    }
1896  }  }
1897    
1898    // For this, we don't care what op you were doing because the answer is now zero
1899    DataAbstract*
1900    DataLazy::zeroedCopy() const
1901    {
1902      return new DataLazy(m_id->zeroedCopy()->getPtr());
1903    }
1904    
1905  // There is no single, natural interpretation of getLength on DataLazy.  // There is no single, natural interpretation of getLength on DataLazy.
1906  // Instances of DataReady can look at the size of their vectors.  // Instances of DataReady can look at the size of their vectors.
# Line 2012  DataLazy::deepCopy() Line 1908  DataLazy::deepCopy()
1908  // 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
1909  // form part of the expression.  // form part of the expression.
1910  // Rather than have people making assumptions, I have disabled the method.  // Rather than have people making assumptions, I have disabled the method.
1911  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1912  DataLazy::getLength() const  DataLazy::getLength() const
1913  {  {
1914    throw DataException("getLength() does not make sense for lazy data.");    throw DataException("getLength() does not make sense for lazy data.");
# Line 2027  DataLazy::getSlice(const DataTypes::Regi Line 1923  DataLazy::getSlice(const DataTypes::Regi
1923    
1924    
1925  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
1926  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1927  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1928                   int dataPointNo)                   int dataPointNo)
1929  {  {
1930    if (m_op==IDENTITY)    if (m_op==IDENTITY)
1931    {    {
1932      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1933    }    }
1934    if (m_readytype!='E')    if (m_readytype!='E')
1935    {    {
1936      collapse();          collapse();
1937      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1938    }    }
1939    // 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
1940    // so we only need to know which child to ask    // so we only need to know which child to ask
1941    if (m_left->m_readytype=='E')    if (m_left->m_readytype=='E')
1942    {    {
1943      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
1944    }    }
1945    else    else
1946    {    {
1947      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
1948    }    }
1949  }  }
1950    
1951  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
1952  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
1953  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1954                   int dataPointNo) const                   int dataPointNo) const
1955  {  {
1956    if (m_op==IDENTITY)    if (m_op==IDENTITY)
1957    {    {
1958      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
1959    }    }
1960    if (m_readytype=='E')    if (m_readytype=='E')
1961    {    {
# Line 2067  DataLazy::getPointOffset(int sampleNo, Line 1963  DataLazy::getPointOffset(int sampleNo,
1963      // so we only need to know which child to ask      // so we only need to know which child to ask
1964      if (m_left->m_readytype=='E')      if (m_left->m_readytype=='E')
1965      {      {
1966      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
1967      }      }
1968      else      else
1969      {      {
1970      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
1971      }      }
1972    }    }
1973    if (m_readytype=='C')    if (m_readytype=='C')
1974    {    {
1975      return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter          return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1976    }    }
1977    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).");
1978  }  }
# Line 2086  DataLazy::getPointOffset(int sampleNo, Line 1982  DataLazy::getPointOffset(int sampleNo,
1982  void  void
1983  DataLazy::setToZero()  DataLazy::setToZero()
1984  {  {
1985  //   DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::RealVectorType v(getNoValues(),0);
1986  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1987  //   m_op=IDENTITY;  //   m_op=IDENTITY;
1988  //   m_right.reset();    //   m_right.reset();  
# Line 2094  DataLazy::setToZero() Line 1990  DataLazy::setToZero()
1990  //   m_readytype='C';  //   m_readytype='C';
1991  //   m_buffsRequired=1;  //   m_buffsRequired=1;
1992    
1993    privdebug=privdebug;  // to stop the compiler complaining about unused privdebug    (void)privdebug;  // to stop the compiler complaining about unused privdebug
1994    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).");
1995  }  }
1996    
1997  bool  bool
1998  DataLazy::actsExpanded() const  DataLazy::actsExpanded() const
1999  {  {
2000      return (m_readytype=='E');          return (m_readytype=='E');
2001  }  }
2002    
2003  }   // end namespace  } // end namespace
2004    

Legend:
Removed from v.4286  
changed lines
  Added in v.6469

  ViewVC Help
Powered by ViewVC 1.1.26