/[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 4264 by jfenwick, Thu Feb 28 11:32:57 2013 UTC trunk/escriptcore/src/DataLazy.cpp revision 5972 by caltinay, Wed Feb 24 04:05:30 2016 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
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
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    #define ESNEEDPYTHON
18    #include "esysUtils/first.h"
19    
20  #include "DataLazy.h"  #include "DataLazy.h"
 #include "esysUtils/Esys_MPI.h"  
 #ifdef _OPENMP  
 #include <omp.h>  
 #endif  
 #include "FunctionSpace.h"  
 #include "DataTypes.h"  
21  #include "Data.h"  #include "Data.h"
22  #include "UnaryFuncs.h"     // for escript::fsign  #include "DataTypes.h"
23    #include "EscriptParams.h"
24    #include "FunctionSpace.h"
25    #include "UnaryFuncs.h"    // for escript::fsign
26  #include "Utils.h"  #include "Utils.h"
27    
28  #include "EscriptParams.h"  #include "esysUtils/Esys_MPI.h"
29    
30    #ifdef _OPENMP
31    #include <omp.h>
32    #endif
33  #ifdef USE_NETCDF  #ifdef USE_NETCDF
34  #include <netcdfcpp.h>  #include <netcdfcpp.h>
35  #endif  #endif
36    
37  #include <iomanip>      // for some fancy formatting in debug  #include <iomanip>              // for some fancy formatting in debug
38    
39    using namespace escript::DataTypes;
40    
41    #define NO_ARG
42    
43  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
44  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
# Line 61  A special operation, IDENTITY, stores an Line 68  A special operation, IDENTITY, stores an
68  This means that all "internal" nodes in the structure are instances of DataLazy.  This means that all "internal" nodes in the structure are instances of DataLazy.
69    
70  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
71  Note that IDENITY is not considered a unary operation.  Note that IDENTITY is not considered a unary operation.
72    
73  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
74  It must however form a DAG (directed acyclic graph).  It must however form a DAG (directed acyclic graph).
# Line 69  I will refer to individual DataLazy obje Line 76  I will refer to individual DataLazy obje
76    
77  Each node also stores:  Each node also stores:
78  - 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
79      evaluated.          evaluated.
80  - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to  - m_buffsrequired ~ the large number of samples which would need to be kept simultaneously in order to
81      evaluate the expression.          evaluate the expression.
82  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
83    
84  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 136  enum ES_opgroup Line 143  enum ES_opgroup
143  {  {
144     G_UNKNOWN,     G_UNKNOWN,
145     G_IDENTITY,     G_IDENTITY,
146     G_BINARY,        // pointwise operations with two arguments     G_BINARY,            // pointwise operations with two arguments
147     G_UNARY,     // pointwise operations with one argument     G_UNARY,             // pointwise operations with one argument
148     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     G_UNARY_P,           // pointwise operations with one argument, requiring a parameter
149     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,            // non-pointwise op with one output
150     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,          // non-pointwise op with one output requiring a parameter
151     G_TENSORPROD,    // general tensor product     G_TENSORPROD,        // general tensor product
152     G_NP1OUT_2P,     // non-pointwise op with one output requiring two params     G_NP1OUT_2P,         // non-pointwise op with one output requiring two params
153     G_REDUCTION,     // non-pointwise unary op with a scalar output     G_REDUCTION,         // non-pointwise unary op with a scalar output
154     G_CONDEVAL     G_CONDEVAL
155  };  };
156    
# Line 151  enum ES_opgroup Line 158  enum ES_opgroup
158    
159    
160  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
161              "sin","cos","tan",                          "sin","cos","tan",
162              "asin","acos","atan","sinh","cosh","tanh","erf",                          "asin","acos","atan","sinh","cosh","tanh","erf",
163              "asinh","acosh","atanh",                          "asinh","acosh","atanh",
164              "log10","log","sign","abs","neg","pos","exp","sqrt",                          "log10","log","sign","abs","neg","pos","exp","sqrt",
165              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",                          "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
166              "symmetric","nonsymmetric",                          "symmetric","nonsymmetric",
167              "prod",                          "prod",
168              "transpose", "trace",                          "transpose", "trace",
169              "swapaxes",                          "swapaxes",
170              "minval", "maxval",                          "minval", "maxval",
171              "condEval"};                          "condEval"};
172  int ES_opcount=44;  int ES_opcount=44;
173  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
174              G_UNARY,G_UNARY,G_UNARY, //10                          G_UNARY,G_UNARY,G_UNARY, //10
175              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17                          G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 17
176              G_UNARY,G_UNARY,G_UNARY,                    // 20                          G_UNARY,G_UNARY,G_UNARY,                                        // 20
177              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,G_UNARY,G_UNARY,        // 28
178              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35                          G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,          // 35
179              G_NP1OUT,G_NP1OUT,                          G_NP1OUT,G_NP1OUT,
180              G_TENSORPROD,                          G_TENSORPROD,
181              G_NP1OUT_P, G_NP1OUT_P,                          G_NP1OUT_P, G_NP1OUT_P,
182              G_NP1OUT_2P,                          G_NP1OUT_2P,
183              G_REDUCTION, G_REDUCTION,                          G_REDUCTION, G_REDUCTION,
184              G_CONDEVAL};                          G_CONDEVAL};
185  inline  inline
186  ES_opgroup  ES_opgroup
187  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 186  getOpgroup(ES_optype op) Line 193  getOpgroup(ES_optype op)
193  FunctionSpace  FunctionSpace
194  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
195  {  {
196      // perhaps this should call interpolate and throw or something?          // perhaps this should call interpolate and throw or something?
197      // maybe we need an interpolate node -          // maybe we need an interpolate node -
198      // 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
199      // programming error exception.          // programming error exception.
200    
201    FunctionSpace l=left->getFunctionSpace();    FunctionSpace l=left->getFunctionSpace();
202    FunctionSpace r=right->getFunctionSpace();    FunctionSpace r=right->getFunctionSpace();
# Line 198  resultFS(DataAbstract_ptr left, DataAbst Line 205  resultFS(DataAbstract_ptr left, DataAbst
205      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
206      if (res==1)      if (res==1)
207      {      {
208      return l;          return l;
209      }      }
210      if (res==-1)      if (res==-1)
211      {      {
212      return r;          return r;
213      }      }
214      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
215    }    }
# Line 214  resultFS(DataAbstract_ptr left, DataAbst Line 221  resultFS(DataAbstract_ptr left, DataAbst
221  DataTypes::ShapeType  DataTypes::ShapeType
222  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
223  {  {
224      if (left->getShape()!=right->getShape())          if (left->getShape()!=right->getShape())
225      {          {
226        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))            if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
227        {            {
228          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.");
229        }            }
230    
231        if (left->getRank()==0)   // we need to allow scalar * anything            if (left->getRank()==0)       // we need to allow scalar * anything
232        {            {
233          return right->getShape();                  return right->getShape();
234        }            }
235        if (right->getRank()==0)            if (right->getRank()==0)
236        {            {
237          return left->getShape();                  return left->getShape();
238        }            }
239        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.");
240      }          }
241      return left->getShape();          return left->getShape();
242  }  }
243    
244  // return the shape for "op left"  // return the shape for "op left"
# Line 239  resultShape(DataAbstract_ptr left, DataA Line 246  resultShape(DataAbstract_ptr left, DataA
246  DataTypes::ShapeType  DataTypes::ShapeType
247  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
248  {  {
249      switch(op)          switch(op)
250      {          {
251          case TRANS:          case TRANS:
252         {            // for the scoping of variables             {                    // for the scoping of variables
253          const DataTypes::ShapeType& s=left->getShape();                  const DataTypes::ShapeType& s=left->getShape();
254          DataTypes::ShapeType sh;                  DataTypes::ShapeType sh;
255          int rank=left->getRank();                  int rank=left->getRank();
256          if (axis_offset<0 || axis_offset>rank)                  if (axis_offset<0 || axis_offset>rank)
257          {                  {
258              stringstream e;              stringstream e;
259              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
260              throw DataException(e.str());              throw DataException(e.str());
261          }          }
262          for (int i=0; i<rank; i++)          for (int i=0; i<rank; i++)
263          {                  {
264             int index = (axis_offset+i)%rank;                     int index = (axis_offset+i)%rank;
265             sh.push_back(s[index]); // Append to new shape             sh.push_back(s[index]); // Append to new shape
266          }          }
267          return sh;                  return sh;
268         }             }
269      break;          break;
270      case TRACE:          case TRACE:
271         {             {
272          int rank=left->getRank();                  int rank=left->getRank();
273          if (rank<2)                  if (rank<2)
274          {                  {
275             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.");
276          }                  }
277          if ((axis_offset>rank-2) || (axis_offset<0))                  if ((axis_offset>rank-2) || (axis_offset<0))
278          {                  {
279             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.");
280          }                  }
281          if (rank==2)                  if (rank==2)
282          {                  {
283             return DataTypes::scalarShape;                     return DataTypes::scalarShape;
284          }                  }
285          else if (rank==3)                  else if (rank==3)
286          {                  {
287             DataTypes::ShapeType sh;                     DataTypes::ShapeType sh;
288                 if (axis_offset==0)                     if (axis_offset==0)
289             {                     {
290                  sh.push_back(left->getShape()[2]);                          sh.push_back(left->getShape()[2]);
291                 }                     }
292                 else     // offset==1                     else         // offset==1
293             {                     {
294              sh.push_back(left->getShape()[0]);                          sh.push_back(left->getShape()[0]);
295                 }                     }
296             return sh;                     return sh;
297          }                  }
298          else if (rank==4)                  else if (rank==4)
299          {                  {
300             DataTypes::ShapeType sh;                     DataTypes::ShapeType sh;
301             const DataTypes::ShapeType& s=left->getShape();                     const DataTypes::ShapeType& s=left->getShape();
302                 if (axis_offset==0)                     if (axis_offset==0)
303             {                     {
304                  sh.push_back(s[2]);                          sh.push_back(s[2]);
305                  sh.push_back(s[3]);                          sh.push_back(s[3]);
306                 }                     }
307                 else if (axis_offset==1)                     else if (axis_offset==1)
308             {                     {
309                  sh.push_back(s[0]);                          sh.push_back(s[0]);
310                  sh.push_back(s[3]);                          sh.push_back(s[3]);
311                 }                     }
312             else     // offset==2                     else         // offset==2
313             {                     {
314              sh.push_back(s[0]);                          sh.push_back(s[0]);
315              sh.push_back(s[1]);                          sh.push_back(s[1]);
316             }                     }
317             return sh;                     return sh;
318          }                  }
319          else        // unknown rank                  else            // unknown rank
320          {                  {
321             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.");
322          }                  }
323         }             }
324      break;          break;
325          default:          default:
326      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)+".");
327      }          }
328  }  }
329    
330  DataTypes::ShapeType  DataTypes::ShapeType
# Line 373  SwapShape(DataAbstract_ptr left, const i Line 380  SwapShape(DataAbstract_ptr left, const i
380  DataTypes::ShapeType  DataTypes::ShapeType
381  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)
382  {  {
383                
384    // Get rank and shape of inputs    // Get rank and shape of inputs
385    int rank0 = left->getRank();    int rank0 = left->getRank();
386    int rank1 = right->getRank();    int rank1 = right->getRank();
# Line 382  GTPShape(DataAbstract_ptr left, DataAbst Line 389  GTPShape(DataAbstract_ptr left, DataAbst
389    
390    // Prepare for the loops of the product and verify compatibility of shapes    // Prepare for the loops of the product and verify compatibility of shapes
391    int start0=0, start1=0;    int start0=0, start1=0;
392    if (transpose == 0)       {}    if (transpose == 0)           {}
393    else if (transpose == 1)  { start0 = axis_offset; }    else if (transpose == 1)      { start0 = axis_offset; }
394    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)      { start1 = rank1-axis_offset; }
395    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"); }
396    
397    if (rank0<axis_offset)    if (rank0<axis_offset)
398    {    {
399      throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");          throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
400    }    }
401    
402    // Adjust the shapes for transpose    // Adjust the shapes for transpose
403    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather    DataTypes::ShapeType tmpShape0(rank0);        // pre-sizing the vectors rather
404    DataTypes::ShapeType tmpShape1(rank1);    // than using push_back    DataTypes::ShapeType tmpShape1(rank1);        // than using push_back
405    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]; }
406    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]; }
407    
408    // Prepare for the loops of the product    // Prepare for the loops of the product
409    SL=1, SM=1, SR=1;    SL=1, SM=1, SR=1;
410    for (int i=0; i<rank0-axis_offset; i++)   {    for (int i=0; i<rank0-axis_offset; i++)       {
411      SL *= tmpShape0[i];      SL *= tmpShape0[i];
412    }    }
413    for (int i=rank0-axis_offset; i<rank0; i++)   {    for (int i=rank0-axis_offset; i<rank0; i++)   {
414      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
415        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
416      }      }
417      SM *= tmpShape0[i];      SM *= tmpShape0[i];
418    }    }
419    for (int i=axis_offset; i<rank1; i++)     {    for (int i=axis_offset; i<rank1; i++)         {
420      SR *= tmpShape1[i];      SR *= tmpShape1[i];
421    }    }
422    
423    // 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)
424    DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      
425    {         // block to limit the scope of out_index    {                     // block to limit the scope of out_index
426       int out_index=0;       int out_index=0;
427       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
428       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 438  GTPShape(DataAbstract_ptr left, DataAbst
438    return shape2;    return shape2;
439  }  }
440    
441  }   // end anonymous namespace  }       // end anonymous namespace
442    
443    
444    
# Line 466  void DataLazy::LazyNodeSetup() Line 473  void DataLazy::LazyNodeSetup()
473    
474  // Creates an identity node  // Creates an identity node
475  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
476      : parent(p->getFunctionSpace(),p->getShape())          : parent(p->getFunctionSpace(),p->getShape())
477      ,m_sampleids(0),          ,m_sampleids(0),
478      m_samples(1)          m_samples(1)
479  {  {
480     if (p->isLazy())     if (p->isLazy())
481     {     {
482      // I don't want identity of Lazy.          // I don't want identity of Lazy.
483      // Question: Why would that be so bad?          // Question: Why would that be so bad?
484      // 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
485      throw DataException("Programmer error - attempt to create identity from a DataLazy.");          throw DataException("Programmer error - attempt to create identity from a DataLazy.");
486     }     }
487     else     else
488     {     {
489      p->makeLazyShared();          p->makeLazyShared();
490      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);          DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
491      makeIdentity(dr);          makeIdentity(dr);
492  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
493     }     }
494  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
495  }  }
496    
497  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
498      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),          : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
499      m_op(op),          m_op(op),
500      m_axis_offset(0),          m_axis_offset(0),
501      m_transpose(0),          m_transpose(0),
502      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
503  {  {
504     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))
505     {     {
506      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.");
507     }     }
508    
509     DataLazy_ptr lleft;     DataLazy_ptr lleft;
510     if (!left->isLazy())     if (!left->isLazy())
511     {     {
512      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
513     }     }
514     else     else
515     {     {
516      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
517     }     }
518     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
519     m_left=lleft;     m_left=lleft;
# Line 520  DataLazy::DataLazy(DataAbstract_ptr left Line 527  DataLazy::DataLazy(DataAbstract_ptr left
527    
528  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
529  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
530      : parent(resultFS(left,right,op), resultShape(left,right,op)),          : parent(resultFS(left,right,op), resultShape(left,right,op)),
531      m_op(op),          m_op(op),
532      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
533  {  {
534  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
535     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
536     {     {
537      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.");
538     }     }
539    
540     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
541     {     {
542      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
543      Data ltemp(left);          Data ltemp(left);
544      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
545      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
546     }     }
547     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
548     {     {
549      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
550      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
551  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
552     }     }
553     left->operandCheck(*right);     left->operandCheck(*right);
554    
555     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
556     {     {
557      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
558  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
559     }     }
560     else     else
561     {     {
562      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
563  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
564     }     }
565     if (right->isLazy())     if (right->isLazy())
566     {     {
567      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
568  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
569     }     }
570     else     else
571     {     {
572      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
573  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
574     }     }
575     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
576     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
577     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
578     {     {
579      m_readytype='E';          m_readytype='E';
580     }     }
581     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
582     {     {
583      m_readytype='T';          m_readytype='T';
584     }     }
585     else     else
586     {     {
587      m_readytype='C';          m_readytype='C';
588     }     }
589     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
590     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 595  LAZYDEBUG(cout << "(3)Lazy created with
595  }  }
596    
597  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)
598      : 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)),
599      m_op(op),          m_op(op),
600      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
601      m_transpose(transpose)          m_transpose(transpose)
602  {  {
603     if ((getOpgroup(op)!=G_TENSORPROD))     if ((getOpgroup(op)!=G_TENSORPROD))
604     {     {
605      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.");
606     }     }
607     if ((transpose>2) || (transpose<0))     if ((transpose>2) || (transpose<0))
608     {     {
609      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");
610     }     }
611     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
612     {     {
613      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
614      Data ltemp(left);          Data ltemp(left);
615      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
616      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
617     }     }
618     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
619     {     {
620      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
621      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
622     }     }
623  //    left->operandCheck(*right);  //    left->operandCheck(*right);
624    
625     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
626     {     {
627      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
628     }     }
629     else     else
630     {     {
631      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
632     }     }
633     if (right->isLazy())     if (right->isLazy())
634     {     {
635      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
636     }     }
637     else     else
638     {     {
639      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
640     }     }
641     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
642     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
643     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
644     {     {
645      m_readytype='E';          m_readytype='E';
646     }     }
647     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
648     {     {
649      m_readytype='T';          m_readytype='T';
650     }     }
651     else     else
652     {     {
653      m_readytype='C';          m_readytype='C';
654     }     }
655     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
656     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 662  LAZYDEBUG(cout << "(4)Lazy created with
662    
663    
664  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
665      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),          : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
666      m_op(op),          m_op(op),
667      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
668      m_transpose(0),          m_transpose(0),
669      m_tol(0)          m_tol(0)
670  {  {
671     if ((getOpgroup(op)!=G_NP1OUT_P))     if ((getOpgroup(op)!=G_NP1OUT_P))
672     {     {
673      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.");
674     }     }
675     DataLazy_ptr lleft;     DataLazy_ptr lleft;
676     if (!left->isLazy())     if (!left->isLazy())
677     {     {
678      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
679     }     }
680     else     else
681     {     {
682      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
683     }     }
684     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
685     m_left=lleft;     m_left=lleft;
# Line 685  LAZYDEBUG(cout << "(5)Lazy created with Line 692  LAZYDEBUG(cout << "(5)Lazy created with
692  }  }
693    
694  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
695      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
696      m_op(op),          m_op(op),
697      m_axis_offset(0),          m_axis_offset(0),
698      m_transpose(0),          m_transpose(0),
699      m_tol(tol)          m_tol(tol)
700  {  {
701     if ((getOpgroup(op)!=G_UNARY_P))     if ((getOpgroup(op)!=G_UNARY_P))
702     {     {
703      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.");
704     }     }
705     DataLazy_ptr lleft;     DataLazy_ptr lleft;
706     if (!left->isLazy())     if (!left->isLazy())
707     {     {
708      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
709     }     }
710     else     else
711     {     {
712      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
713     }     }
714     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
715     m_left=lleft;     m_left=lleft;
# Line 716  LAZYDEBUG(cout << "(6)Lazy created with Line 723  LAZYDEBUG(cout << "(6)Lazy created with
723    
724    
725  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)
726      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),          : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
727      m_op(op),          m_op(op),
728      m_axis_offset(axis0),          m_axis_offset(axis0),
729      m_transpose(axis1),          m_transpose(axis1),
730      m_tol(0)          m_tol(0)
731  {  {
732     if ((getOpgroup(op)!=G_NP1OUT_2P))     if ((getOpgroup(op)!=G_NP1OUT_2P))
733     {     {
734      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.");
735     }     }
736     DataLazy_ptr lleft;     DataLazy_ptr lleft;
737     if (!left->isLazy())     if (!left->isLazy())
738     {     {
739      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
740     }     }
741     else     else
742     {     {
743      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
744     }     }
745     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
746     m_left=lleft;     m_left=lleft;
# Line 751  namespace Line 758  namespace
758    
759      inline int max3(int a, int b, int c)      inline int max3(int a, int b, int c)
760      {      {
761      int t=(a>b?a:b);          int t=(a>b?a:b);
762      return (t>c?t:c);          return (t>c?t:c);
763    
764      }      }
765  }  }
766    
767  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*/)
768      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
769      m_op(CONDEVAL),          m_op(CONDEVAL),
770      m_axis_offset(0),          m_axis_offset(0),
771      m_transpose(0),          m_transpose(0),
772      m_tol(0)          m_tol(0)
773  {  {
774    
775     DataLazy_ptr lmask;     DataLazy_ptr lmask;
# Line 770  DataLazy::DataLazy(DataAbstract_ptr mask Line 777  DataLazy::DataLazy(DataAbstract_ptr mask
777     DataLazy_ptr lright;     DataLazy_ptr lright;
778     if (!mask->isLazy())     if (!mask->isLazy())
779     {     {
780      lmask=DataLazy_ptr(new DataLazy(mask));          lmask=DataLazy_ptr(new DataLazy(mask));
781     }     }
782     else     else
783     {     {
784      lmask=dynamic_pointer_cast<DataLazy>(mask);          lmask=dynamic_pointer_cast<DataLazy>(mask);
785     }     }
786     if (!left->isLazy())     if (!left->isLazy())
787     {     {
788      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
789     }     }
790     else     else
791     {     {
792      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
793     }     }
794     if (!right->isLazy())     if (!right->isLazy())
795     {     {
796      lright=DataLazy_ptr(new DataLazy(right));          lright=DataLazy_ptr(new DataLazy(right));
797     }     }
798     else     else
799     {     {
800      lright=dynamic_pointer_cast<DataLazy>(right);          lright=dynamic_pointer_cast<DataLazy>(right);
801     }     }
802     m_readytype=lmask->m_readytype;     m_readytype=lmask->m_readytype;
803     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))
804     {     {
805      throw DataException("Programmer Error - condEval arguments must have the same readytype");          throw DataException("Programmer Error - condEval arguments must have the same readytype");
806     }     }
807     m_left=lleft;     m_left=lleft;
808     m_right=lright;     m_right=lright;
# Line 822  DataLazy::~DataLazy() Line 829  DataLazy::~DataLazy()
829    For reasons of efficiency do not call this method on DataExpanded nodes.    For reasons of efficiency do not call this method on DataExpanded nodes.
830  */  */
831  DataReady_ptr  DataReady_ptr
832  DataLazy::collapseToReady()  DataLazy::collapseToReady() const
833  {  {
834    if (m_readytype=='E')    if (m_readytype=='E')
835    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
836      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
837    }    }
838    if (m_op==IDENTITY)    if (m_op==IDENTITY)
# Line 843  DataLazy::collapseToReady() Line 850  DataLazy::collapseToReady()
850    switch(m_op)    switch(m_op)
851    {    {
852      case ADD:      case ADD:
853      result=left+right;          result=left+right;
854      break;          break;
855      case SUB:            case SUB:          
856      result=left-right;          result=left-right;
857      break;          break;
858      case MUL:            case MUL:          
859      result=left*right;          result=left*right;
860      break;          break;
861      case DIV:            case DIV:          
862      result=left/right;          result=left/right;
863      break;          break;
864      case SIN:      case SIN:
865      result=left.sin();            result=left.sin();      
866      break;          break;
867      case COS:      case COS:
868      result=left.cos();          result=left.cos();
869      break;          break;
870      case TAN:      case TAN:
871      result=left.tan();          result=left.tan();
872      break;          break;
873      case ASIN:      case ASIN:
874      result=left.asin();          result=left.asin();
875      break;          break;
876      case ACOS:      case ACOS:
877      result=left.acos();          result=left.acos();
878      break;          break;
879      case ATAN:      case ATAN:
880      result=left.atan();          result=left.atan();
881      break;          break;
882      case SINH:      case SINH:
883      result=left.sinh();          result=left.sinh();
884      break;          break;
885      case COSH:      case COSH:
886      result=left.cosh();          result=left.cosh();
887      break;          break;
888      case TANH:      case TANH:
889      result=left.tanh();          result=left.tanh();
890      break;          break;
891      case ERF:      case ERF:
892      result=left.erf();          result=left.erf();
893      break;          break;
894     case ASINH:     case ASINH:
895      result=left.asinh();          result=left.asinh();
896      break;          break;
897     case ACOSH:     case ACOSH:
898      result=left.acosh();          result=left.acosh();
899      break;          break;
900     case ATANH:     case ATANH:
901      result=left.atanh();          result=left.atanh();
902      break;          break;
903      case LOG10:      case LOG10:
904      result=left.log10();          result=left.log10();
905      break;          break;
906      case LOG:      case LOG:
907      result=left.log();          result=left.log();
908      break;          break;
909      case SIGN:      case SIGN:
910      result=left.sign();          result=left.sign();
911      break;          break;
912      case ABS:      case ABS:
913      result=left.abs();          result=left.abs();
914      break;          break;
915      case NEG:      case NEG:
916      result=left.neg();          result=left.neg();
917      break;          break;
918      case POS:      case POS:
919      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
920      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
921      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
922      break;          break;
923      case EXP:      case EXP:
924      result=left.exp();          result=left.exp();
925      break;          break;
926      case SQRT:      case SQRT:
927      result=left.sqrt();          result=left.sqrt();
928      break;          break;
929      case RECIP:      case RECIP:
930      result=left.oneOver();          result=left.oneOver();
931      break;          break;
932      case GZ:      case GZ:
933      result=left.wherePositive();          result=left.wherePositive();
934      break;          break;
935      case LZ:      case LZ:
936      result=left.whereNegative();          result=left.whereNegative();
937      break;          break;
938      case GEZ:      case GEZ:
939      result=left.whereNonNegative();          result=left.whereNonNegative();
940      break;          break;
941      case LEZ:      case LEZ:
942      result=left.whereNonPositive();          result=left.whereNonPositive();
943      break;          break;
944      case NEZ:      case NEZ:
945      result=left.whereNonZero(m_tol);          result=left.whereNonZero(m_tol);
946      break;          break;
947      case EZ:      case EZ:
948      result=left.whereZero(m_tol);          result=left.whereZero(m_tol);
949      break;          break;
950      case SYM:      case SYM:
951      result=left.symmetric();          result=left.symmetric();
952      break;          break;
953      case NSYM:      case NSYM:
954      result=left.nonsymmetric();          result=left.nonsymmetric();
955      break;          break;
956      case PROD:      case PROD:
957      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);          result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
958      break;          break;
959      case TRANS:      case TRANS:
960      result=left.transpose(m_axis_offset);          result=left.transpose(m_axis_offset);
961      break;          break;
962      case TRACE:      case TRACE:
963      result=left.trace(m_axis_offset);          result=left.trace(m_axis_offset);
964      break;          break;
965      case SWAP:      case SWAP:
966      result=left.swapaxes(m_axis_offset, m_transpose);          result=left.swapaxes(m_axis_offset, m_transpose);
967      break;          break;
968      case MINVAL:      case MINVAL:
969      result=left.minval();          result=left.minval();
970      break;          break;
971      case MAXVAL:      case MAXVAL:
972      result=left.minval();          result=left.minval();
973      break;          break;
974      default:      default:
975      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)+".");
976    }    }
977    return result.borrowReadyPtr();    return result.borrowReadyPtr();
978  }  }
# Line 977  DataLazy::collapseToReady() Line 984  DataLazy::collapseToReady()
984     the purpose of using DataLazy in the first place).     the purpose of using DataLazy in the first place).
985  */  */
986  void  void
987  DataLazy::collapse()  DataLazy::collapse() const
988  {  {
989    if (m_op==IDENTITY)    if (m_op==IDENTITY)
990    {    {
991      return;          return;
992    }    }
993    if (m_readytype=='E')    if (m_readytype=='E')
994    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
995      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
996    }    }
997    m_id=collapseToReady();    m_id=collapseToReady();
# Line 997  DataLazy::collapse() Line 1004  DataLazy::collapse()
1004    
1005    
1006  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1007      for (int j=0;j<onumsteps;++j)\          for (int j=0;j<onumsteps;++j)\
1008      {\          {\
1009        for (int i=0;i<numsteps;++i,resultp+=resultStep) \            for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1010        { \            { \
1011  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1012  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1013           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \               tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1014  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1015           lroffset+=leftstep; \               lroffset+=leftstep; \
1016           rroffset+=rightstep; \               rroffset+=rightstep; \
1017        }\            }\
1018        lroffset+=oleftstep;\            lroffset+=oleftstep;\
1019        rroffset+=orightstep;\            rroffset+=orightstep;\
1020      }          }
1021    
1022    
1023  // The result will be stored in m_samples  // The result will be stored in m_samples
1024  // 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
1025  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1026  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
1027  {  {
1028  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1029      // 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
1030    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1031    {    {
1032      collapse();          collapse();
1033    }    }
1034    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
1035    {    {
1036      const ValueType& vec=m_id->getVectorRO();      const ValueType& vec=m_id->getVectorRO();
1037      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
# Line 1043  if (&x<stackend[omp_get_thread_num()]) Line 1050  if (&x<stackend[omp_get_thread_num()])
1050    }    }
1051    if (m_sampleids[tid]==sampleNo)    if (m_sampleids[tid]==sampleNo)
1052    {    {
1053      roffset=tid*m_samplesize;          roffset=tid*m_samplesize;
1054      return &(m_samples);        // sample is already resolved          return &(m_samples);            // sample is already resolved
1055    }    }
1056    m_sampleids[tid]=sampleNo;    m_sampleids[tid]=sampleNo;
1057    
# Line 1064  if (&x<stackend[omp_get_thread_num()]) Line 1071  if (&x<stackend[omp_get_thread_num()])
1071    }    }
1072  }  }
1073    
1074  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1075  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
1076  {  {
1077      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1078      // 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
1079      // processing single points.          // processing single points.
1080      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1081    if (m_readytype!='E')    if (m_readytype!='E')
1082    {    {
1083      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 1086  DataLazy::resolveNodeUnary(int tid, int
1086    {    {
1087      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1088    }    }
1089    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1090    const double* left=&((*leftres)[roffset]);    const double* left=&((*leftres)[roffset]);
1091    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1092    double* result=&(m_samples[roffset]);    double* result=&(m_samples[roffset]);
1093    switch (m_op)    switch (m_op)
1094    {    {
1095      case SIN:        case SIN:  
1096      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
1097      break;          break;
1098      case COS:      case COS:
1099      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
1100      break;          break;
1101      case TAN:      case TAN:
1102      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
1103      break;          break;
1104      case ASIN:      case ASIN:
1105      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1106      break;          break;
1107      case ACOS:      case ACOS:
1108      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1109      break;          break;
1110      case ATAN:      case ATAN:
1111      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1112      break;          break;
1113      case SINH:      case SINH:
1114      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1115      break;          break;
1116      case COSH:      case COSH:
1117      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1118      break;          break;
1119      case TANH:      case TANH:
1120      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1121      break;          break;
1122      case ERF:      case ERF:
1123  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1124      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");          throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1125  #else  #else
1126      tensor_unary_operation(m_samplesize, left, result, ::erf);          tensor_unary_operation(m_samplesize, left, result, ::erf);
1127      break;          break;
1128  #endif  #endif
1129     case ASINH:     case ASINH:
1130  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1131      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);          tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1132  #else  #else
1133      tensor_unary_operation(m_samplesize, left, result, ::asinh);          tensor_unary_operation(m_samplesize, left, result, ::asinh);
1134  #endif    #endif  
1135      break;          break;
1136     case ACOSH:     case ACOSH:
1137  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1138      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);          tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1139  #else  #else
1140      tensor_unary_operation(m_samplesize, left, result, ::acosh);          tensor_unary_operation(m_samplesize, left, result, ::acosh);
1141  #endif    #endif  
1142      break;          break;
1143     case ATANH:     case ATANH:
1144  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1145      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);          tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1146  #else  #else
1147      tensor_unary_operation(m_samplesize, left, result, ::atanh);          tensor_unary_operation(m_samplesize, left, result, ::atanh);
1148  #endif    #endif  
1149      break;          break;
1150      case LOG10:      case LOG10:
1151      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1152      break;          break;
1153      case LOG:      case LOG:
1154      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1155      break;          break;
1156      case SIGN:      case SIGN:
1157      tensor_unary_operation(m_samplesize, left, result, escript::fsign);          tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1158      break;          break;
1159      case ABS:      case ABS:
1160      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1161      break;          break;
1162      case NEG:      case NEG:
1163      tensor_unary_operation(m_samplesize, left, result, negate<double>());          tensor_unary_operation(m_samplesize, left, result, negate<double>());
1164      break;          break;
1165      case POS:      case POS:
1166      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
1167      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
1168      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
1169      break;          break;
1170      case EXP:      case EXP:
1171      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1172      break;          break;
1173      case SQRT:      case SQRT:
1174      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1175      break;          break;
1176      case RECIP:      case RECIP:
1177      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));          tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1178      break;          break;
1179      case GZ:      case GZ:
1180      tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));          tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1181      break;          break;
1182      case LZ:      case LZ:
1183      tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));          tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1184      break;          break;
1185      case GEZ:      case GEZ:
1186      tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));          tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1187      break;          break;
1188      case LEZ:      case LEZ:
1189      tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));          tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1190      break;          break;
1191  // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently  // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1192      case NEZ:      case NEZ:
1193      tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));          tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1194      break;          break;
1195      case EZ:      case EZ:
1196      tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));          tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1197      break;          break;
1198    
1199      default:      default:
1200      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1201    }    }
1202    return &(m_samples);    return &(m_samples);
1203  }  }
1204    
1205    
1206  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1207  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1208  {  {
1209      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1210      // 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
1211      // processing single points.          // processing single points.
1212      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1213    if (m_readytype!='E')    if (m_readytype!='E')
1214    {    {
1215      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 1219  DataLazy::resolveNodeReduction(int tid,
1219      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1220    }    }
1221    size_t loffset=0;    size_t loffset=0;
1222    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1223    
1224    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1225    unsigned int ndpps=getNumDPPSample();    unsigned int ndpps=getNumDPPSample();
# Line 1221  DataLazy::resolveNodeReduction(int tid, Line 1228  DataLazy::resolveNodeReduction(int tid,
1228    switch (m_op)    switch (m_op)
1229    {    {
1230      case MINVAL:      case MINVAL:
1231      {          {
1232        for (unsigned int z=0;z<ndpps;++z)            for (unsigned int z=0;z<ndpps;++z)
1233        {            {
1234          FMin op;              FMin op;
1235          *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());              *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1236          loffset+=psize;              loffset+=psize;
1237          result++;              result++;
1238        }            }
1239      }          }
1240      break;          break;
1241      case MAXVAL:      case MAXVAL:
1242      {          {
1243        for (unsigned int z=0;z<ndpps;++z)            for (unsigned int z=0;z<ndpps;++z)
1244        {            {
1245        FMax op;            FMax op;
1246        *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);            *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1247        loffset+=psize;            loffset+=psize;
1248        result++;            result++;
1249        }            }
1250      }          }
1251      break;          break;
1252      default:      default:
1253      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1254    }    }
1255    return &(m_samples);    return &(m_samples);
1256  }  }
1257    
1258  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1259  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1260  {  {
1261      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1262      // 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
1263      // processing single points.          // processing single points.
1264    if (m_readytype!='E')    if (m_readytype!='E')
1265    {    {
1266      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 1272  DataLazy::resolveNodeNP1OUT(int tid, int Line 1279  DataLazy::resolveNodeNP1OUT(int tid, int
1279    switch (m_op)    switch (m_op)
1280    {    {
1281      case SYM:      case SYM:
1282      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1283      {          {
1284          DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1285          subroffset+=step;              subroffset+=step;
1286          offset+=step;              offset+=step;
1287      }          }
1288      break;          break;
1289      case NSYM:      case NSYM:
1290      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1291      {          {
1292          DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1293          subroffset+=step;              subroffset+=step;
1294          offset+=step;              offset+=step;
1295      }          }
1296      break;          break;
1297      default:      default:
1298      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1299    }    }
1300    return &m_samples;    return &m_samples;
1301  }  }
1302    
1303  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1304  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1305  {  {
1306      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1307      // 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
1308      // processing single points.          // processing single points.
1309    if (m_readytype!='E')    if (m_readytype!='E')
1310    {    {
1311      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 1319  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1326  DataLazy::resolveNodeNP1OUT_P(int tid, i
1326    switch (m_op)    switch (m_op)
1327    {    {
1328      case TRACE:      case TRACE:
1329      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1330      {          {
1331              DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);              DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1332          subroffset+=instep;              subroffset+=instep;
1333          offset+=outstep;              offset+=outstep;
1334      }          }
1335      break;          break;
1336      case TRANS:      case TRANS:
1337      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1338      {          {
1339              DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);              DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1340          subroffset+=instep;              subroffset+=instep;
1341          offset+=outstep;              offset+=outstep;
1342      }          }
1343      break;          break;
1344      default:      default:
1345      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1346    }    }
1347    return &m_samples;    return &m_samples;
1348  }  }
1349    
1350    
1351  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1352  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1353  {  {
1354    if (m_readytype!='E')    if (m_readytype!='E')
1355    {    {
# Line 1364  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1371  DataLazy::resolveNodeNP1OUT_2P(int tid,
1371    switch (m_op)    switch (m_op)
1372    {    {
1373      case SWAP:      case SWAP:
1374      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1375      {          {
1376              DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);              DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1377          subroffset+=instep;              subroffset+=instep;
1378          offset+=outstep;              offset+=outstep;
1379      }          }
1380      break;          break;
1381      default:      default:
1382      throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1383    }    }
1384    return &m_samples;    return &m_samples;
1385  }  }
1386    
1387  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1388  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1389  {  {
1390    if (m_readytype!='E')    if (m_readytype!='E')
1391    {    {
# Line 1394  DataLazy::resolveNodeCondEval(int tid, i Line 1401  DataLazy::resolveNodeCondEval(int tid, i
1401    const ValueType* srcres=0;    const ValueType* srcres=0;
1402    if ((*maskres)[subroffset]>0)    if ((*maskres)[subroffset]>0)
1403    {    {
1404      srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1405    }    }
1406    else    else
1407    {    {
1408      srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1409    }    }
1410    
1411    // Now we need to copy the result    // Now we need to copy the result
# Line 1406  DataLazy::resolveNodeCondEval(int tid, i Line 1413  DataLazy::resolveNodeCondEval(int tid, i
1413    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1414    for (int i=0;i<m_samplesize;++i)    for (int i=0;i<m_samplesize;++i)
1415    {    {
1416      m_samples[roffset+i]=(*srcres)[subroffset+i];            m_samples[roffset+i]=(*srcres)[subroffset+i];  
1417    }    }
1418    
1419    return &m_samples;    return &m_samples;
# Line 1421  DataLazy::resolveNodeCondEval(int tid, i Line 1428  DataLazy::resolveNodeCondEval(int tid, i
1428  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1429  // For example, 2+Vector.  // For example, 2+Vector.
1430  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1431  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1432  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1433  {  {
1434  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1435    
1436    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
1437      // first work out which of the children are expanded          // first work out which of the children are expanded
1438    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1439    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1440    if (!leftExp && !rightExp)    if (!leftExp && !rightExp)
1441    {    {
1442      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'.");
1443    }    }
1444    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1445    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1446    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1447    {    {
1448      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.");
1449    }    }
1450    size_t leftsize=m_left->getNoValues();    size_t leftsize=m_left->getNoValues();
1451    size_t rightsize=m_right->getNoValues();    size_t rightsize=m_right->getNoValues();
1452    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
1453    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
1454    int rightstep=0;    int rightstep=0;
1455    int numsteps=0;       // total number of steps for the inner loop    int numsteps=0;               // total number of steps for the inner loop
1456    int oleftstep=0;  // the o variables refer to the outer loop    int oleftstep=0;      // the o variables refer to the outer loop
1457    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
1458    int onumsteps=1;    int onumsteps=1;
1459        
1460    bool LES=(leftExp && leftScalar); // Left is an expanded scalar    bool LES=(leftExp && leftScalar);     // Left is an expanded scalar
1461    bool RES=(rightExp && rightScalar);    bool RES=(rightExp && rightScalar);
1462    bool LS=(!leftExp && leftScalar); // left is a single scalar    bool LS=(!leftExp && leftScalar);     // left is a single scalar
1463    bool RS=(!rightExp && rightScalar);    bool RS=(!rightExp && rightScalar);
1464    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1465    bool RN=(!rightExp && !rightScalar);    bool RN=(!rightExp && !rightScalar);
1466    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1467    bool REN=(rightExp && !rightScalar);    bool REN=(rightExp && !rightScalar);
1468    
1469    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
1470    {    {
1471      chunksize=m_left->getNumDPPSample()*leftsize;          chunksize=m_left->getNumDPPSample()*leftsize;
1472      leftstep=0;          leftstep=0;
1473      rightstep=0;          rightstep=0;
1474      numsteps=1;          numsteps=1;
1475    }    }
1476    else if (LES || RES)    else if (LES || RES)
1477    {    {
1478      chunksize=1;          chunksize=1;
1479      if (LES)        // left is an expanded scalar          if (LES)                // left is an expanded scalar
1480      {          {
1481          if (RS)                  if (RS)
1482          {                  {
1483             leftstep=1;                     leftstep=1;
1484             rightstep=0;                     rightstep=0;
1485             numsteps=m_left->getNumDPPSample();                     numsteps=m_left->getNumDPPSample();
1486          }                  }
1487          else        // RN or REN                  else            // RN or REN
1488          {                  {
1489             leftstep=0;                     leftstep=0;
1490             oleftstep=1;                     oleftstep=1;
1491             rightstep=1;                     rightstep=1;
1492             orightstep=(RN ? -(int)rightsize : 0);                     orightstep=(RN ? -(int)rightsize : 0);
1493             numsteps=rightsize;                     numsteps=rightsize;
1494             onumsteps=m_left->getNumDPPSample();                     onumsteps=m_left->getNumDPPSample();
1495          }                  }
1496      }          }
1497      else        // right is an expanded scalar          else            // right is an expanded scalar
1498      {          {
1499          if (LS)                  if (LS)
1500          {                  {
1501             rightstep=1;                     rightstep=1;
1502             leftstep=0;                     leftstep=0;
1503             numsteps=m_right->getNumDPPSample();                     numsteps=m_right->getNumDPPSample();
1504          }                  }
1505          else                  else
1506          {                  {
1507             rightstep=0;                     rightstep=0;
1508             orightstep=1;                     orightstep=1;
1509             leftstep=1;                     leftstep=1;
1510             oleftstep=(LN ? -(int)leftsize : 0);                     oleftstep=(LN ? -(int)leftsize : 0);
1511             numsteps=leftsize;                     numsteps=leftsize;
1512             onumsteps=m_right->getNumDPPSample();                     onumsteps=m_right->getNumDPPSample();
1513          }                  }
1514      }          }
1515    }    }
1516    else  // this leaves (LEN, RS), (LEN, RN) and their transposes    else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1517    {    {
1518      if (LEN)    // and Right will be a single value          if (LEN)        // and Right will be a single value
1519      {          {
1520          chunksize=rightsize;                  chunksize=rightsize;
1521          leftstep=rightsize;                  leftstep=rightsize;
1522          rightstep=0;                  rightstep=0;
1523          numsteps=m_left->getNumDPPSample();                  numsteps=m_left->getNumDPPSample();
1524          if (RS)                  if (RS)
1525          {                  {
1526             numsteps*=leftsize;                     numsteps*=leftsize;
1527          }                  }
1528      }          }
1529      else    // REN          else    // REN
1530      {          {
1531          chunksize=leftsize;                  chunksize=leftsize;
1532          rightstep=leftsize;                  rightstep=leftsize;
1533          leftstep=0;                  leftstep=0;
1534          numsteps=m_right->getNumDPPSample();                  numsteps=m_right->getNumDPPSample();
1535          if (LS)                  if (LS)
1536          {                  {
1537             numsteps*=rightsize;                     numsteps*=rightsize;
1538          }                  }
1539      }          }
1540    }    }
1541    
1542    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1543      // Get the values of sub-expressions          // Get the values of sub-expressions
1544    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1545    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1546  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1547  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;)
# Line 1549  LAZYDEBUG(cout << "Right res["<< rroffse Line 1556  LAZYDEBUG(cout << "Right res["<< rroffse
1556    
1557    
1558    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1559    double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved    double* resultp=&(m_samples[roffset]);                // results are stored at the vector offset we received
1560    switch(m_op)    switch(m_op)
1561    {    {
1562      case ADD:      case ADD:
1563          PROC_OP(NO_ARG,plus<double>());          PROC_OP(NO_ARG,plus<double>());
1564      break;          break;
1565      case SUB:      case SUB:
1566      PROC_OP(NO_ARG,minus<double>());          PROC_OP(NO_ARG,minus<double>());
1567      break;          break;
1568      case MUL:      case MUL:
1569      PROC_OP(NO_ARG,multiplies<double>());          PROC_OP(NO_ARG,multiplies<double>());
1570      break;          break;
1571      case DIV:      case DIV:
1572      PROC_OP(NO_ARG,divides<double>());          PROC_OP(NO_ARG,divides<double>());
1573      break;          break;
1574      case POW:      case POW:
1575         PROC_OP(double (double,double),::pow);         PROC_OP(double (double,double),::pow);
1576      break;          break;
1577      default:      default:
1578      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1579    }    }
1580  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1581    return &m_samples;    return &m_samples;
# Line 1578  LAZYDEBUG(cout << "Result res[" << roffs Line 1585  LAZYDEBUG(cout << "Result res[" << roffs
1585  // 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
1586  // 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.
1587  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1588  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1589  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1590  {  {
1591  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1592    
1593    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
1594      // first work out which of the children are expanded          // first work out which of the children are expanded
1595    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1596    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1597    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1598    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
1599    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1600    
1601    int resultStep=getNoValues();    int resultStep=getNoValues();
# Line 1611  LAZYDEBUG(cout << "m_samplesize=" << m_s Line 1618  LAZYDEBUG(cout << "m_samplesize=" << m_s
1618  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1619  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1620    
1621    double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(m_samples[offset]);         // results are stored at the vector offset we received
1622    switch(m_op)    switch(m_op)
1623    {    {
1624      case PROD:      case PROD:
1625      for (int i=0;i<steps;++i,resultp+=resultStep)          for (int i=0;i<steps;++i,resultp+=resultStep)
1626      {          {
1627            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1628            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1629    
1630  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1631  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1632    
1633            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);
1634    
1635        lroffset+=leftStep;            lroffset+=leftStep;
1636        rroffset+=rightStep;            rroffset+=rightStep;
1637      }          }
1638      break;          break;
1639      default:      default:
1640      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1641    }    }
1642    roffset=offset;    roffset=offset;
1643    return &m_samples;    return &m_samples;
1644  }  }
1645    
1646    
1647  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1648  DataLazy::resolveSample(int sampleNo, size_t& roffset)  DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1649  {  {
1650  #ifdef _OPENMP  #ifdef _OPENMP
1651      int tid=omp_get_thread_num();          int tid=omp_get_thread_num();
1652  #else  #else
1653      int tid=0;          int tid=0;
1654  #endif  #endif
1655    
1656  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1657      stackstart[tid]=&tid;          stackstart[tid]=&tid;
1658      stackend[tid]=&tid;          stackend[tid]=&tid;
1659      const DataTypes::ValueType* r=resolveNodeSample(tid, sampleNo, roffset);          const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1660      size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];          size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1661      #pragma omp critical          #pragma omp critical
1662      if (d>maxstackuse)          if (d>maxstackuse)
1663      {          {
1664  cout << "Max resolve Stack use " << d << endl;  cout << "Max resolve Stack use " << d << endl;
1665          maxstackuse=d;                  maxstackuse=d;
1666      }          }
1667      return r;          return r;
1668  #else  #else
1669      return resolveNodeSample(tid, sampleNo, roffset);          return resolveNodeSample(tid, sampleNo, roffset);
1670  #endif  #endif
1671  }  }
1672    
# Line 1669  void Line 1676  void
1676  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
1677  {  {
1678     if (m_op==IDENTITY)     if (m_op==IDENTITY)
1679      return;          return;
1680     DataReady_ptr p=resolveNodeWorker();     DataReady_ptr p=resolveNodeWorker();
1681     makeIdentity(p);     makeIdentity(p);
1682  }  }
# Line 1706  DataLazy::resolveGroupWorker(std::vector Line 1713  DataLazy::resolveGroupWorker(std::vector
1713  {  {
1714    if (dats.empty())    if (dats.empty())
1715    {    {
1716      return;          return;
1717    }    }
1718    vector<DataLazy*> work;    vector<DataLazy*> work;
1719    FunctionSpace fs=dats[0]->getFunctionSpace();    FunctionSpace fs=dats[0]->getFunctionSpace();
1720    bool match=true;    bool match=true;
1721    for (int i=dats.size()-1;i>=0;--i)    for (int i=dats.size()-1;i>=0;--i)
1722    {    {
1723      if (dats[i]->m_readytype!='E')          if (dats[i]->m_readytype!='E')
1724      {          {
1725          dats[i]->collapse();                  dats[i]->collapse();
1726      }          }
1727      if (dats[i]->m_op!=IDENTITY)          if (dats[i]->m_op!=IDENTITY)
1728      {          {
1729          work.push_back(dats[i]);                  work.push_back(dats[i]);
1730          if (fs!=dats[i]->getFunctionSpace())                  if (fs!=dats[i]->getFunctionSpace())
1731          {                  {
1732              match=false;                          match=false;
1733          }                  }
1734      }          }
1735    }    }
1736    if (work.empty())    if (work.empty())
1737    {    {
1738      return;     // no work to do          return;         // no work to do
1739    }    }
1740    if (match)    // all functionspaces match.  Yes I realise this is overly strict    if (match)    // all functionspaces match.  Yes I realise this is overly strict
1741    {     // 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
1742          // all the other functionspaces match.                  // all the other functionspaces match.
1743      vector<DataExpanded*> dep;          vector<DataExpanded*> dep;
1744      vector<ValueType*> vecs;          vector<ValueType*> vecs;
1745      for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1746      {          {
1747          dep.push_back(new DataExpanded(fs,work[i]->getShape(), ValueType(work[i]->getNoValues())));                  dep.push_back(new DataExpanded(fs,work[i]->getShape(), ValueType(work[i]->getNoValues())));
1748          vecs.push_back(&(dep[i]->getVectorRW()));                  vecs.push_back(&(dep[i]->getVectorRW()));
1749      }          }
1750      int totalsamples=work[0]->getNumSamples();          int totalsamples=work[0]->getNumSamples();
1751      const ValueType* res=0; // Storage for answer          const ValueType* res=0; // Storage for answer
1752      int sample;          int sample;
1753      #pragma omp parallel private(sample, res)          #pragma omp parallel private(sample, res)
1754      {          {
1755          size_t roffset=0;              size_t roffset=0;
1756          #pragma omp for schedule(static)              #pragma omp for schedule(static)
1757          for (sample=0;sample<totalsamples;++sample)              for (sample=0;sample<totalsamples;++sample)
1758          {              {
1759          roffset=0;                  roffset=0;
1760          int j;                  int j;
1761          for (j=work.size()-1;j>=0;--j)                  for (j=work.size()-1;j>=0;--j)
1762          {                  {
1763  #ifdef _OPENMP  #ifdef _OPENMP
1764                  res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);                      res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1765  #else  #else
1766                  res=work[j]->resolveNodeSample(0,sample,roffset);                      res=work[j]->resolveNodeSample(0,sample,roffset);
1767  #endif  #endif
1768                  DataVector::size_type outoffset=dep[j]->getPointOffset(sample,0);                      RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1769                  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));
1770          }                  }
1771          }              }
1772      }          }
1773      // 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
1774      for (int i=work.size()-1;i>=0;--i)          for (int i=work.size()-1;i>=0;--i)
1775      {          {
1776          work[i]->makeIdentity(boost::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));              work[i]->makeIdentity(boost::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1777      }          }
1778    }    }
1779    else  // functionspaces do not match    else  // functionspaces do not match
1780    {    {
1781      for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1782      {          {
1783          work[i]->resolveToIdentity();                  work[i]->resolveToIdentity();
1784      }          }
1785    }    }
1786  }  }
1787    
# Line 1784  DataLazy::resolveGroupWorker(std::vector Line 1791  DataLazy::resolveGroupWorker(std::vector
1791  DataReady_ptr  DataReady_ptr
1792  DataLazy::resolveNodeWorker()  DataLazy::resolveNodeWorker()
1793  {  {
1794    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
1795    {    {
1796      collapse();      collapse();
1797    }    }
1798    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.
1799    {    {
1800      return m_id;      return m_id;
1801    }    }
1802      // 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'
1803    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1804    ValueType& resvec=result->getVectorRW();    ValueType& resvec=result->getVectorRW();
1805    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1806    
1807    int sample;    int sample;
1808    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1809    const ValueType* res=0;   // Storage for answer    const ValueType* res=0;       // Storage for answer
1810  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1811    #pragma omp parallel private(sample,res)    #pragma omp parallel private(sample,res)
1812    {    {
1813      size_t roffset=0;          size_t roffset=0;
1814  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1815      stackstart[omp_get_thread_num()]=&roffset;          stackstart[omp_get_thread_num()]=&roffset;
1816      stackend[omp_get_thread_num()]=&roffset;          stackend[omp_get_thread_num()]=&roffset;
1817  #endif  #endif
1818      #pragma omp for schedule(static)          #pragma omp for schedule(static)
1819      for (sample=0;sample<totalsamples;++sample)          for (sample=0;sample<totalsamples;++sample)
1820      {          {
1821          roffset=0;                  roffset=0;
1822  #ifdef _OPENMP  #ifdef _OPENMP
1823              res=resolveNodeSample(omp_get_thread_num(),sample,roffset);                  res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1824  #else  #else
1825              res=resolveNodeSample(0,sample,roffset);                  res=resolveNodeSample(0,sample,roffset);
1826  #endif  #endif
1827  LAZYDEBUG(cout << "Sample #" << sample << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1828  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1829              DataVector::size_type outoffset=result->getPointOffset(sample,0);                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1830              memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));                  memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1831      }          }
1832    }    }
1833  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1834    for (int i=0;i<getNumberOfThreads();++i)    for (int i=0;i<getNumberOfThreads();++i)
1835    {    {
1836      size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);          size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1837  //  cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;  //      cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1838      if (r>maxstackuse)          if (r>maxstackuse)
1839      {          {
1840          maxstackuse=r;                  maxstackuse=r;
1841      }          }
1842    }    }
1843    cout << "Max resolve Stack use=" << maxstackuse << endl;    cout << "Max resolve Stack use=" << maxstackuse << endl;
1844  #endif  #endif
# Line 1845  DataLazy::toString() const Line 1852  DataLazy::toString() const
1852    oss << "Lazy Data: [depth=" << m_height<< "] ";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1853    switch (escriptParams.getLAZY_STR_FMT())    switch (escriptParams.getLAZY_STR_FMT())
1854    {    {
1855    case 1:   // tree format    case 1:       // tree format
1856      oss << endl;          oss << endl;
1857      intoTreeString(oss,"");          intoTreeString(oss,"");
1858      break;          break;
1859    case 2:   // just the depth    case 2:       // just the depth
1860      break;          break;
1861    default:    default:
1862      intoString(oss);          intoString(oss);
1863      break;          break;
1864    }    }
1865    return oss.str();    return oss.str();
1866  }  }
# Line 1866  DataLazy::intoString(ostringstream& oss) Line 1873  DataLazy::intoString(ostringstream& oss)
1873    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1874    {    {
1875    case G_IDENTITY:    case G_IDENTITY:
1876      if (m_id->isExpanded())          if (m_id->isExpanded())
1877      {          {
1878         oss << "E";             oss << "E";
1879      }          }
1880      else if (m_id->isTagged())          else if (m_id->isTagged())
1881      {          {
1882        oss << "T";            oss << "T";
1883      }          }
1884      else if (m_id->isConstant())          else if (m_id->isConstant())
1885      {          {
1886        oss << "C";            oss << "C";
1887      }          }
1888      else          else
1889      {          {
1890        oss << "?";            oss << "?";
1891      }          }
1892      oss << '@' << m_id.get();          oss << '@' << m_id.get();
1893      break;          break;
1894    case G_BINARY:    case G_BINARY:
1895      oss << '(';          oss << '(';
1896      m_left->intoString(oss);          m_left->intoString(oss);
1897      oss << ' ' << opToString(m_op) << ' ';          oss << ' ' << opToString(m_op) << ' ';
1898      m_right->intoString(oss);          m_right->intoString(oss);
1899      oss << ')';          oss << ')';
1900      break;          break;
1901    case G_UNARY:    case G_UNARY:
1902    case G_UNARY_P:    case G_UNARY_P:
1903    case G_NP1OUT:    case G_NP1OUT:
1904    case G_NP1OUT_P:    case G_NP1OUT_P:
1905    case G_REDUCTION:    case G_REDUCTION:
1906      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1907      m_left->intoString(oss);          m_left->intoString(oss);
1908      oss << ')';          oss << ')';
1909      break;          break;
1910    case G_TENSORPROD:    case G_TENSORPROD:
1911      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1912      m_left->intoString(oss);          m_left->intoString(oss);
1913      oss << ", ";          oss << ", ";
1914      m_right->intoString(oss);          m_right->intoString(oss);
1915      oss << ')';          oss << ')';
1916      break;          break;
1917    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1918      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1919      m_left->intoString(oss);          m_left->intoString(oss);
1920      oss << ", " << m_axis_offset << ", " << m_transpose;          oss << ", " << m_axis_offset << ", " << m_transpose;
1921      oss << ')';          oss << ')';
1922      break;          break;
1923    case G_CONDEVAL:    case G_CONDEVAL:
1924      oss << opToString(m_op)<< '(' ;          oss << opToString(m_op)<< '(' ;
1925      m_mask->intoString(oss);          m_mask->intoString(oss);
1926      oss << " ? ";          oss << " ? ";
1927      m_left->intoString(oss);          m_left->intoString(oss);
1928      oss << " : ";          oss << " : ";
1929      m_right->intoString(oss);          m_right->intoString(oss);
1930      oss << ')';          oss << ')';
1931      break;          break;
1932    default:    default:
1933      oss << "UNKNOWN";          oss << "UNKNOWN";
1934    }    }
1935  }  }
1936    
# Line 1935  DataLazy::intoTreeString(ostringstream& Line 1942  DataLazy::intoTreeString(ostringstream&
1942    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1943    {    {
1944    case G_IDENTITY:    case G_IDENTITY:
1945      if (m_id->isExpanded())          if (m_id->isExpanded())
1946      {          {
1947         oss << "E";             oss << "E";
1948      }          }
1949      else if (m_id->isTagged())          else if (m_id->isTagged())
1950      {          {
1951        oss << "T";            oss << "T";
1952      }          }
1953      else if (m_id->isConstant())          else if (m_id->isConstant())
1954      {          {
1955        oss << "C";            oss << "C";
1956      }          }
1957      else          else
1958      {          {
1959        oss << "?";            oss << "?";
1960      }          }
1961      oss << '@' << m_id.get() << endl;          oss << '@' << m_id.get() << endl;
1962      break;          break;
1963    case G_BINARY:    case G_BINARY:
1964      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1965      indent+='.';          indent+='.';
1966      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1967      m_right->intoTreeString(oss, indent);          m_right->intoTreeString(oss, indent);
1968      break;          break;
1969    case G_UNARY:    case G_UNARY:
1970    case G_UNARY_P:    case G_UNARY_P:
1971    case G_NP1OUT:    case G_NP1OUT:
1972    case G_NP1OUT_P:    case G_NP1OUT_P:
1973    case G_REDUCTION:    case G_REDUCTION:
1974      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1975      indent+='.';          indent+='.';
1976      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1977      break;          break;
1978    case G_TENSORPROD:    case G_TENSORPROD:
1979      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1980      indent+='.';          indent+='.';
1981      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1982      m_right->intoTreeString(oss, indent);          m_right->intoTreeString(oss, indent);
1983      break;          break;
1984    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1985      oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;          oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1986      indent+='.';          indent+='.';
1987      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1988      break;          break;
1989    default:    default:
1990      oss << "UNKNOWN";          oss << "UNKNOWN";
1991    }    }
1992  }  }
1993    
1994    
1995  DataAbstract*  DataAbstract*
1996  DataLazy::deepCopy()  DataLazy::deepCopy() const
1997  {  {
1998    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1999    {    {
2000    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
2001    case G_UNARY:    case G_UNARY:
2002    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
2003    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);
2004    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);
2005    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);
2006    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);
2007    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);
2008    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);
2009    default:    default:
2010      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)+".");
2011    }    }
2012  }  }
2013    
# Line 2012  DataLazy::deepCopy() Line 2019  DataLazy::deepCopy()
2019  // 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
2020  // form part of the expression.  // form part of the expression.
2021  // Rather than have people making assumptions, I have disabled the method.  // Rather than have people making assumptions, I have disabled the method.
2022  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
2023  DataLazy::getLength() const  DataLazy::getLength() const
2024  {  {
2025    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 2034  DataLazy::getSlice(const DataTypes::Regi
2034    
2035    
2036  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
2037  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
2038  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
2039                   int dataPointNo)                   int dataPointNo)
2040  {  {
2041    if (m_op==IDENTITY)    if (m_op==IDENTITY)
2042    {    {
2043      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
2044    }    }
2045    if (m_readytype!='E')    if (m_readytype!='E')
2046    {    {
2047      collapse();          collapse();
2048      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
2049    }    }
2050    // 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
2051    // so we only need to know which child to ask    // so we only need to know which child to ask
2052    if (m_left->m_readytype=='E')    if (m_left->m_readytype=='E')
2053    {    {
2054      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
2055    }    }
2056    else    else
2057    {    {
2058      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
2059    }    }
2060  }  }
2061    
2062  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
2063  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
2064  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
2065                   int dataPointNo) const                   int dataPointNo) const
2066  {  {
2067    if (m_op==IDENTITY)    if (m_op==IDENTITY)
2068    {    {
2069      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
2070    }    }
2071    if (m_readytype=='E')    if (m_readytype=='E')
2072    {    {
# Line 2067  DataLazy::getPointOffset(int sampleNo, Line 2074  DataLazy::getPointOffset(int sampleNo,
2074      // so we only need to know which child to ask      // so we only need to know which child to ask
2075      if (m_left->m_readytype=='E')      if (m_left->m_readytype=='E')
2076      {      {
2077      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
2078      }      }
2079      else      else
2080      {      {
2081      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
2082      }      }
2083    }    }
2084    if (m_readytype=='C')    if (m_readytype=='C')
2085    {    {
2086      return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter          return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
2087    }    }
2088    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).");
2089  }  }
# Line 2086  DataLazy::getPointOffset(int sampleNo, Line 2093  DataLazy::getPointOffset(int sampleNo,
2093  void  void
2094  DataLazy::setToZero()  DataLazy::setToZero()
2095  {  {
2096  //   DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::RealVectorType v(getNoValues(),0);
2097  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2098  //   m_op=IDENTITY;  //   m_op=IDENTITY;
2099  //   m_right.reset();    //   m_right.reset();  
# Line 2094  DataLazy::setToZero() Line 2101  DataLazy::setToZero()
2101  //   m_readytype='C';  //   m_readytype='C';
2102  //   m_buffsRequired=1;  //   m_buffsRequired=1;
2103    
2104    privdebug=privdebug;  // to stop the compiler complaining about unused privdebug    (void)privdebug;  // to stop the compiler complaining about unused privdebug
2105    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).");
2106  }  }
2107    
2108  bool  bool
2109  DataLazy::actsExpanded() const  DataLazy::actsExpanded() const
2110  {  {
2111      return (m_readytype=='E');          return (m_readytype=='E');
2112  }  }
2113    
2114  }   // end namespace  }       // end namespace
2115    

Legend:
Removed from v.4264  
changed lines
  Added in v.5972

  ViewVC Help
Powered by ViewVC 1.1.26