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

revision 5938 by jfenwick, Thu Feb 18 06:30:35 2016 UTC revision 5972 by caltinay, Wed Feb 24 04:05:30 2016 UTC
# Line 18  Line 18 
18  #include "esysUtils/first.h"  #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;  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)
45  namespace  namespace
# Line 74  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 large number of samples which would need to be kept simultaneously in order to  - m_buffsrequired ~ the large number of samples which would need to be kept simultaneously in order to
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 141  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 156  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 191  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 203  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 219  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 244  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 378  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 387  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 436  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 471  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 525  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 593  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 660  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 690  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 721  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 756  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 775  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 830  DataReady_ptr Line 832  DataReady_ptr
832  DataLazy::collapseToReady() const  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 848  DataLazy::collapseToReady() const Line 850  DataLazy::collapseToReady() const
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 986  DataLazy::collapse() const Line 988  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 1002  DataLazy::collapse() const Line 1004  DataLazy::collapse() const
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
# Line 1024  const DataTypes::RealVectorType* Line 1026  const DataTypes::RealVectorType*
1026  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const  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 1048  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 1072  if (&x<stackend[omp_get_thread_num()]) Line 1074  if (&x<stackend[omp_get_thread_num()])
1074  const DataTypes::RealVectorType*  const DataTypes::RealVectorType*
1075  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const  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 1090  DataLazy::resolveNodeUnary(int tid, int Line 1092  DataLazy::resolveNodeUnary(int tid, int
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  }  }
# Line 1204  DataLazy::resolveNodeUnary(int tid, int Line 1206  DataLazy::resolveNodeUnary(int tid, int
1206  const DataTypes::RealVectorType*  const DataTypes::RealVectorType*
1207  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const  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 1226  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  }  }
# Line 1256  DataLazy::resolveNodeReduction(int tid, Line 1258  DataLazy::resolveNodeReduction(int tid,
1258  const DataTypes::RealVectorType*  const DataTypes::RealVectorType*
1259  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const  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 1277  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  }  }
# Line 1301  DataLazy::resolveNodeNP1OUT(int tid, int Line 1303  DataLazy::resolveNodeNP1OUT(int tid, int
1303  const DataTypes::RealVectorType*  const DataTypes::RealVectorType*
1304  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const  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 1324  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  }  }
# Line 1369  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  }  }
# Line 1399  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 1411  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 1431  DataLazy::resolveNodeBinary(int tid, int Line 1433  DataLazy::resolveNodeBinary(int tid, int
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 1554  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 received    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 1588  DataLazy::resolveNodeTProd(int tid, int Line 1590  DataLazy::resolveNodeTProd(int tid, int
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 1616  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 received    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;
# Line 1646  const DataTypes::RealVectorType* Line 1648  const DataTypes::RealVectorType*
1648  DataLazy::resolveSample(int sampleNo, size_t& roffset) const  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::RealVectorType* 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 1674  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 1711  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                  RealVectorType::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(RealVectorType::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 1789  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              RealVectorType::size_type outoffset=result->getPointOffset(sample,0);                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1830              memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::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 1850  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 1871  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 1940  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    
# Line 1996  DataLazy::deepCopy() const Line 1998  DataLazy::deepCopy() const
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 2038  DataLazy::getPointOffset(int sampleNo, Line 2040  DataLazy::getPointOffset(int sampleNo,
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    
# Line 2064  DataLazy::getPointOffset(int sampleNo, Line 2066  DataLazy::getPointOffset(int sampleNo,
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 2072  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 2106  DataLazy::setToZero() Line 2108  DataLazy::setToZero()
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.5938  
changed lines
  Added in v.5972

  ViewVC Help
Powered by ViewVC 1.1.26