/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/clazy/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/escript/src/DataLazy.cpp revision 4264 by jfenwick, Thu Feb 28 11:32:57 2013 UTC trunk/escriptcore/src/DataLazy.cpp revision 6001 by caltinay, Tue Mar 1 05:01:49 2016 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2016 by The University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
   
17  #include "DataLazy.h"  #include "DataLazy.h"
 #include "esysUtils/Esys_MPI.h"  
 #ifdef _OPENMP  
 #include <omp.h>  
 #endif  
 #include "FunctionSpace.h"  
 #include "DataTypes.h"  
18  #include "Data.h"  #include "Data.h"
19  #include "UnaryFuncs.h"     // for escript::fsign  #include "DataTypes.h"
 #include "Utils.h"  
   
20  #include "EscriptParams.h"  #include "EscriptParams.h"
21    #include "FunctionSpace.h"
22    #include "UnaryFuncs.h"    // for escript::fsign
23    #include "Utils.h"
24    
25  #ifdef USE_NETCDF  #ifdef USE_NETCDF
26  #include <netcdfcpp.h>  #include <netcdfcpp.h>
27  #endif  #endif
28    
29  #include <iomanip>      // for some fancy formatting in debug  #include <iomanip> // for some fancy formatting in debug
30    
31    using namespace escript::DataTypes;
32    
33    #define NO_ARG
34    
35  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
36  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
# Line 61  A special operation, IDENTITY, stores an Line 60  A special operation, IDENTITY, stores an
60  This means that all "internal" nodes in the structure are instances of DataLazy.  This means that all "internal" nodes in the structure are instances of DataLazy.
61    
62  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
63  Note that IDENITY is not considered a unary operation.  Note that IDENTITY is not considered a unary operation.
64    
65  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
66  It must however form a DAG (directed acyclic graph).  It must however form a DAG (directed acyclic graph).
# Line 69  I will refer to individual DataLazy obje Line 68  I will refer to individual DataLazy obje
68    
69  Each node also stores:  Each node also stores:
70  - 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
71      evaluated.          evaluated.
72  - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to  - m_buffsrequired ~ the large number of samples which would need to be kept simultaneously in order to
73      evaluate the expression.          evaluate the expression.
74  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
75    
76  When a new node is created, the above values are computed based on the values in the child nodes.  When a new node is created, the above values are computed based on the values in the child nodes.
# Line 136  enum ES_opgroup Line 135  enum ES_opgroup
135  {  {
136     G_UNKNOWN,     G_UNKNOWN,
137     G_IDENTITY,     G_IDENTITY,
138     G_BINARY,        // pointwise operations with two arguments     G_BINARY,            // pointwise operations with two arguments
139     G_UNARY,     // pointwise operations with one argument     G_UNARY,             // pointwise operations with one argument
140     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     G_UNARY_P,           // pointwise operations with one argument, requiring a parameter
141     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,            // non-pointwise op with one output
142     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,          // non-pointwise op with one output requiring a parameter
143     G_TENSORPROD,    // general tensor product     G_TENSORPROD,        // general tensor product
144     G_NP1OUT_2P,     // non-pointwise op with one output requiring two params     G_NP1OUT_2P,         // non-pointwise op with one output requiring two params
145     G_REDUCTION,     // non-pointwise unary op with a scalar output     G_REDUCTION,         // non-pointwise unary op with a scalar output
146     G_CONDEVAL     G_CONDEVAL
147  };  };
148    
# Line 151  enum ES_opgroup Line 150  enum ES_opgroup
150    
151    
152  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
153              "sin","cos","tan",                          "sin","cos","tan",
154              "asin","acos","atan","sinh","cosh","tanh","erf",                          "asin","acos","atan","sinh","cosh","tanh","erf",
155              "asinh","acosh","atanh",                          "asinh","acosh","atanh",
156              "log10","log","sign","abs","neg","pos","exp","sqrt",                          "log10","log","sign","abs","neg","pos","exp","sqrt",
157              "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",
158              "symmetric","nonsymmetric",                          "symmetric","nonsymmetric",
159              "prod",                          "prod",
160              "transpose", "trace",                          "transpose", "trace",
161              "swapaxes",                          "swapaxes",
162              "minval", "maxval",                          "minval", "maxval",
163              "condEval"};                          "condEval"};
164  int ES_opcount=44;  int ES_opcount=44;
165  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,
166              G_UNARY,G_UNARY,G_UNARY, //10                          G_UNARY,G_UNARY,G_UNARY, //10
167              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
168              G_UNARY,G_UNARY,G_UNARY,                    // 20                          G_UNARY,G_UNARY,G_UNARY,                                        // 20
169              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
170              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
171              G_NP1OUT,G_NP1OUT,                          G_NP1OUT,G_NP1OUT,
172              G_TENSORPROD,                          G_TENSORPROD,
173              G_NP1OUT_P, G_NP1OUT_P,                          G_NP1OUT_P, G_NP1OUT_P,
174              G_NP1OUT_2P,                          G_NP1OUT_2P,
175              G_REDUCTION, G_REDUCTION,                          G_REDUCTION, G_REDUCTION,
176              G_CONDEVAL};                          G_CONDEVAL};
177  inline  inline
178  ES_opgroup  ES_opgroup
179  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 186  getOpgroup(ES_optype op) Line 185  getOpgroup(ES_optype op)
185  FunctionSpace  FunctionSpace
186  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
187  {  {
188      // perhaps this should call interpolate and throw or something?          // perhaps this should call interpolate and throw or something?
189      // maybe we need an interpolate node -          // maybe we need an interpolate node -
190      // 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
191      // programming error exception.          // programming error exception.
192    
193    FunctionSpace l=left->getFunctionSpace();    FunctionSpace l=left->getFunctionSpace();
194    FunctionSpace r=right->getFunctionSpace();    FunctionSpace r=right->getFunctionSpace();
# Line 198  resultFS(DataAbstract_ptr left, DataAbst Line 197  resultFS(DataAbstract_ptr left, DataAbst
197      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
198      if (res==1)      if (res==1)
199      {      {
200      return l;          return l;
201      }      }
202      if (res==-1)      if (res==-1)
203      {      {
204      return r;          return r;
205      }      }
206      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");      throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
207    }    }
# Line 214  resultFS(DataAbstract_ptr left, DataAbst Line 213  resultFS(DataAbstract_ptr left, DataAbst
213  DataTypes::ShapeType  DataTypes::ShapeType
214  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
215  {  {
216      if (left->getShape()!=right->getShape())          if (left->getShape()!=right->getShape())
217      {          {
218        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))            if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
219        {            {
220          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.");
221        }            }
222    
223        if (left->getRank()==0)   // we need to allow scalar * anything            if (left->getRank()==0)       // we need to allow scalar * anything
224        {            {
225          return right->getShape();                  return right->getShape();
226        }            }
227        if (right->getRank()==0)            if (right->getRank()==0)
228        {            {
229          return left->getShape();                  return left->getShape();
230        }            }
231        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.");
232      }          }
233      return left->getShape();          return left->getShape();
234  }  }
235    
236  // return the shape for "op left"  // return the shape for "op left"
# Line 239  resultShape(DataAbstract_ptr left, DataA Line 238  resultShape(DataAbstract_ptr left, DataA
238  DataTypes::ShapeType  DataTypes::ShapeType
239  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
240  {  {
241      switch(op)          switch(op)
242      {          {
243          case TRANS:          case TRANS:
244         {            // for the scoping of variables             {                    // for the scoping of variables
245          const DataTypes::ShapeType& s=left->getShape();                  const DataTypes::ShapeType& s=left->getShape();
246          DataTypes::ShapeType sh;                  DataTypes::ShapeType sh;
247          int rank=left->getRank();                  int rank=left->getRank();
248          if (axis_offset<0 || axis_offset>rank)                  if (axis_offset<0 || axis_offset>rank)
249          {                  {
250              stringstream e;              stringstream e;
251              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
252              throw DataException(e.str());              throw DataException(e.str());
253          }          }
254          for (int i=0; i<rank; i++)          for (int i=0; i<rank; i++)
255          {                  {
256             int index = (axis_offset+i)%rank;                     int index = (axis_offset+i)%rank;
257             sh.push_back(s[index]); // Append to new shape             sh.push_back(s[index]); // Append to new shape
258          }          }
259          return sh;                  return sh;
260         }             }
261      break;          break;
262      case TRACE:          case TRACE:
263         {             {
264          int rank=left->getRank();                  int rank=left->getRank();
265          if (rank<2)                  if (rank<2)
266          {                  {
267             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.");
268          }                  }
269          if ((axis_offset>rank-2) || (axis_offset<0))                  if ((axis_offset>rank-2) || (axis_offset<0))
270          {                  {
271             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.");
272          }                  }
273          if (rank==2)                  if (rank==2)
274          {                  {
275             return DataTypes::scalarShape;                     return DataTypes::scalarShape;
276          }                  }
277          else if (rank==3)                  else if (rank==3)
278          {                  {
279             DataTypes::ShapeType sh;                     DataTypes::ShapeType sh;
280                 if (axis_offset==0)                     if (axis_offset==0)
281             {                     {
282                  sh.push_back(left->getShape()[2]);                          sh.push_back(left->getShape()[2]);
283                 }                     }
284                 else     // offset==1                     else         // offset==1
285             {                     {
286              sh.push_back(left->getShape()[0]);                          sh.push_back(left->getShape()[0]);
287                 }                     }
288             return sh;                     return sh;
289          }                  }
290          else if (rank==4)                  else if (rank==4)
291          {                  {
292             DataTypes::ShapeType sh;                     DataTypes::ShapeType sh;
293             const DataTypes::ShapeType& s=left->getShape();                     const DataTypes::ShapeType& s=left->getShape();
294                 if (axis_offset==0)                     if (axis_offset==0)
295             {                     {
296                  sh.push_back(s[2]);                          sh.push_back(s[2]);
297                  sh.push_back(s[3]);                          sh.push_back(s[3]);
298                 }                     }
299                 else if (axis_offset==1)                     else if (axis_offset==1)
300             {                     {
301                  sh.push_back(s[0]);                          sh.push_back(s[0]);
302                  sh.push_back(s[3]);                          sh.push_back(s[3]);
303                 }                     }
304             else     // offset==2                     else         // offset==2
305             {                     {
306              sh.push_back(s[0]);                          sh.push_back(s[0]);
307              sh.push_back(s[1]);                          sh.push_back(s[1]);
308             }                     }
309             return sh;                     return sh;
310          }                  }
311          else        // unknown rank                  else            // unknown rank
312          {                  {
313             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.");
314          }                  }
315         }             }
316      break;          break;
317          default:          default:
318      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)+".");
319      }          }
320  }  }
321    
322  DataTypes::ShapeType  DataTypes::ShapeType
# Line 373  SwapShape(DataAbstract_ptr left, const i Line 372  SwapShape(DataAbstract_ptr left, const i
372  DataTypes::ShapeType  DataTypes::ShapeType
373  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)
374  {  {
375                
376    // Get rank and shape of inputs    // Get rank and shape of inputs
377    int rank0 = left->getRank();    int rank0 = left->getRank();
378    int rank1 = right->getRank();    int rank1 = right->getRank();
# Line 382  GTPShape(DataAbstract_ptr left, DataAbst Line 381  GTPShape(DataAbstract_ptr left, DataAbst
381    
382    // Prepare for the loops of the product and verify compatibility of shapes    // Prepare for the loops of the product and verify compatibility of shapes
383    int start0=0, start1=0;    int start0=0, start1=0;
384    if (transpose == 0)       {}    if (transpose == 0)           {}
385    else if (transpose == 1)  { start0 = axis_offset; }    else if (transpose == 1)      { start0 = axis_offset; }
386    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)      { start1 = rank1-axis_offset; }
387    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"); }
388    
389    if (rank0<axis_offset)    if (rank0<axis_offset)
390    {    {
391      throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");          throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
392    }    }
393    
394    // Adjust the shapes for transpose    // Adjust the shapes for transpose
395    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather    DataTypes::ShapeType tmpShape0(rank0);        // pre-sizing the vectors rather
396    DataTypes::ShapeType tmpShape1(rank1);    // than using push_back    DataTypes::ShapeType tmpShape1(rank1);        // than using push_back
397    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]; }
398    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]; }
399    
400    // Prepare for the loops of the product    // Prepare for the loops of the product
401    SL=1, SM=1, SR=1;    SL=1, SM=1, SR=1;
402    for (int i=0; i<rank0-axis_offset; i++)   {    for (int i=0; i<rank0-axis_offset; i++)       {
403      SL *= tmpShape0[i];      SL *= tmpShape0[i];
404    }    }
405    for (int i=rank0-axis_offset; i<rank0; i++)   {    for (int i=rank0-axis_offset; i<rank0; i++)   {
406      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {      if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
407        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");        throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
408      }      }
409      SM *= tmpShape0[i];      SM *= tmpShape0[i];
410    }    }
411    for (int i=axis_offset; i<rank1; i++)     {    for (int i=axis_offset; i<rank1; i++)         {
412      SR *= tmpShape1[i];      SR *= tmpShape1[i];
413    }    }
414    
415    // 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)
416    DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      
417    {         // block to limit the scope of out_index    {                     // block to limit the scope of out_index
418       int out_index=0;       int out_index=0;
419       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
420       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
# Line 431  GTPShape(DataAbstract_ptr left, DataAbst Line 430  GTPShape(DataAbstract_ptr left, DataAbst
430    return shape2;    return shape2;
431  }  }
432    
433  }   // end anonymous namespace  }       // end anonymous namespace
434    
435    
436    
# Line 466  void DataLazy::LazyNodeSetup() Line 465  void DataLazy::LazyNodeSetup()
465    
466  // Creates an identity node  // Creates an identity node
467  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
468      : parent(p->getFunctionSpace(),p->getShape())          : parent(p->getFunctionSpace(),p->getShape())
469      ,m_sampleids(0),          ,m_sampleids(0),
470      m_samples(1)          m_samples(1)
471  {  {
472     if (p->isLazy())     if (p->isLazy())
473     {     {
474      // I don't want identity of Lazy.          // I don't want identity of Lazy.
475      // Question: Why would that be so bad?          // Question: Why would that be so bad?
476      // 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
477      throw DataException("Programmer error - attempt to create identity from a DataLazy.");          throw DataException("Programmer error - attempt to create identity from a DataLazy.");
478     }     }
479     else     else
480     {     {
481      p->makeLazyShared();          p->makeLazyShared();
482      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);          DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
483      makeIdentity(dr);          makeIdentity(dr);
484  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
485     }     }
486  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
487  }  }
488    
489  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
490      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),          : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
491      m_op(op),          m_op(op),
492      m_axis_offset(0),          m_axis_offset(0),
493      m_transpose(0),          m_transpose(0),
494      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
495  {  {
496     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))
497     {     {
498      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.");
499     }     }
500    
501     DataLazy_ptr lleft;     DataLazy_ptr lleft;
502     if (!left->isLazy())     if (!left->isLazy())
503     {     {
504      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
505     }     }
506     else     else
507     {     {
508      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
509     }     }
510     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
511     m_left=lleft;     m_left=lleft;
# Line 520  DataLazy::DataLazy(DataAbstract_ptr left Line 519  DataLazy::DataLazy(DataAbstract_ptr left
519    
520  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
521  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
522      : parent(resultFS(left,right,op), resultShape(left,right,op)),          : parent(resultFS(left,right,op), resultShape(left,right,op)),
523      m_op(op),          m_op(op),
524      m_SL(0), m_SM(0), m_SR(0)          m_SL(0), m_SM(0), m_SR(0)
525  {  {
526  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
527     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
528     {     {
529      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.");
530     }     }
531    
532     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
533     {     {
534      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
535      Data ltemp(left);          Data ltemp(left);
536      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
537      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
538     }     }
539     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
540     {     {
541      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
542      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
543  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
544     }     }
545     left->operandCheck(*right);     left->operandCheck(*right);
546    
547     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
548     {     {
549      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
550  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
551     }     }
552     else     else
553     {     {
554      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
555  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
556     }     }
557     if (right->isLazy())     if (right->isLazy())
558     {     {
559      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
560  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
561     }     }
562     else     else
563     {     {
564      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
565  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
566     }     }
567     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
568     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
569     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
570     {     {
571      m_readytype='E';          m_readytype='E';
572     }     }
573     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
574     {     {
575      m_readytype='T';          m_readytype='T';
576     }     }
577     else     else
578     {     {
579      m_readytype='C';          m_readytype='C';
580     }     }
581     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
582     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
# Line 588  LAZYDEBUG(cout << "(3)Lazy created with Line 587  LAZYDEBUG(cout << "(3)Lazy created with
587  }  }
588    
589  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)
590      : 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)),
591      m_op(op),          m_op(op),
592      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
593      m_transpose(transpose)          m_transpose(transpose)
594  {  {
595     if ((getOpgroup(op)!=G_TENSORPROD))     if ((getOpgroup(op)!=G_TENSORPROD))
596     {     {
597      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.");
598     }     }
599     if ((transpose>2) || (transpose<0))     if ((transpose>2) || (transpose<0))
600     {     {
601      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");
602     }     }
603     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
604     {     {
605      FunctionSpace fs=getFunctionSpace();          FunctionSpace fs=getFunctionSpace();
606      Data ltemp(left);          Data ltemp(left);
607      Data tmp(ltemp,fs);          Data tmp(ltemp,fs);
608      left=tmp.borrowDataPtr();          left=tmp.borrowDataPtr();
609     }     }
610     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
611     {     {
612      Data tmp(Data(right),getFunctionSpace());          Data tmp(Data(right),getFunctionSpace());
613      right=tmp.borrowDataPtr();          right=tmp.borrowDataPtr();
614     }     }
615  //    left->operandCheck(*right);  //    left->operandCheck(*right);
616    
617     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
618     {     {
619      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
620     }     }
621     else     else
622     {     {
623      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
624     }     }
625     if (right->isLazy())     if (right->isLazy())
626     {     {
627      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
628     }     }
629     else     else
630     {     {
631      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
632     }     }
633     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
634     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
635     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
636     {     {
637      m_readytype='E';          m_readytype='E';
638     }     }
639     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
640     {     {
641      m_readytype='T';          m_readytype='T';
642     }     }
643     else     else
644     {     {
645      m_readytype='C';          m_readytype='C';
646     }     }
647     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
648     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
# Line 655  LAZYDEBUG(cout << "(4)Lazy created with Line 654  LAZYDEBUG(cout << "(4)Lazy created with
654    
655    
656  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
657      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),          : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
658      m_op(op),          m_op(op),
659      m_axis_offset(axis_offset),          m_axis_offset(axis_offset),
660      m_transpose(0),          m_transpose(0),
661      m_tol(0)          m_tol(0)
662  {  {
663     if ((getOpgroup(op)!=G_NP1OUT_P))     if ((getOpgroup(op)!=G_NP1OUT_P))
664     {     {
665      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.");
666     }     }
667     DataLazy_ptr lleft;     DataLazy_ptr lleft;
668     if (!left->isLazy())     if (!left->isLazy())
669     {     {
670      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
671     }     }
672     else     else
673     {     {
674      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
675     }     }
676     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
677     m_left=lleft;     m_left=lleft;
# Line 685  LAZYDEBUG(cout << "(5)Lazy created with Line 684  LAZYDEBUG(cout << "(5)Lazy created with
684  }  }
685    
686  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
687      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
688      m_op(op),          m_op(op),
689      m_axis_offset(0),          m_axis_offset(0),
690      m_transpose(0),          m_transpose(0),
691      m_tol(tol)          m_tol(tol)
692  {  {
693     if ((getOpgroup(op)!=G_UNARY_P))     if ((getOpgroup(op)!=G_UNARY_P))
694     {     {
695      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.");
696     }     }
697     DataLazy_ptr lleft;     DataLazy_ptr lleft;
698     if (!left->isLazy())     if (!left->isLazy())
699     {     {
700      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
701     }     }
702     else     else
703     {     {
704      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
705     }     }
706     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
707     m_left=lleft;     m_left=lleft;
# Line 716  LAZYDEBUG(cout << "(6)Lazy created with Line 715  LAZYDEBUG(cout << "(6)Lazy created with
715    
716    
717  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)
718      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),          : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
719      m_op(op),          m_op(op),
720      m_axis_offset(axis0),          m_axis_offset(axis0),
721      m_transpose(axis1),          m_transpose(axis1),
722      m_tol(0)          m_tol(0)
723  {  {
724     if ((getOpgroup(op)!=G_NP1OUT_2P))     if ((getOpgroup(op)!=G_NP1OUT_2P))
725     {     {
726      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.");
727     }     }
728     DataLazy_ptr lleft;     DataLazy_ptr lleft;
729     if (!left->isLazy())     if (!left->isLazy())
730     {     {
731      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
732     }     }
733     else     else
734     {     {
735      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
736     }     }
737     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
738     m_left=lleft;     m_left=lleft;
# Line 751  namespace Line 750  namespace
750    
751      inline int max3(int a, int b, int c)      inline int max3(int a, int b, int c)
752      {      {
753      int t=(a>b?a:b);          int t=(a>b?a:b);
754      return (t>c?t:c);          return (t>c?t:c);
755    
756      }      }
757  }  }
758    
759  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*/)
760      : parent(left->getFunctionSpace(), left->getShape()),          : parent(left->getFunctionSpace(), left->getShape()),
761      m_op(CONDEVAL),          m_op(CONDEVAL),
762      m_axis_offset(0),          m_axis_offset(0),
763      m_transpose(0),          m_transpose(0),
764      m_tol(0)          m_tol(0)
765  {  {
766    
767     DataLazy_ptr lmask;     DataLazy_ptr lmask;
# Line 770  DataLazy::DataLazy(DataAbstract_ptr mask Line 769  DataLazy::DataLazy(DataAbstract_ptr mask
769     DataLazy_ptr lright;     DataLazy_ptr lright;
770     if (!mask->isLazy())     if (!mask->isLazy())
771     {     {
772      lmask=DataLazy_ptr(new DataLazy(mask));          lmask=DataLazy_ptr(new DataLazy(mask));
773     }     }
774     else     else
775     {     {
776      lmask=dynamic_pointer_cast<DataLazy>(mask);          lmask=dynamic_pointer_cast<DataLazy>(mask);
777     }     }
778     if (!left->isLazy())     if (!left->isLazy())
779     {     {
780      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
781     }     }
782     else     else
783     {     {
784      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
785     }     }
786     if (!right->isLazy())     if (!right->isLazy())
787     {     {
788      lright=DataLazy_ptr(new DataLazy(right));          lright=DataLazy_ptr(new DataLazy(right));
789     }     }
790     else     else
791     {     {
792      lright=dynamic_pointer_cast<DataLazy>(right);          lright=dynamic_pointer_cast<DataLazy>(right);
793     }     }
794     m_readytype=lmask->m_readytype;     m_readytype=lmask->m_readytype;
795     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))
796     {     {
797      throw DataException("Programmer Error - condEval arguments must have the same readytype");          throw DataException("Programmer Error - condEval arguments must have the same readytype");
798     }     }
799     m_left=lleft;     m_left=lleft;
800     m_right=lright;     m_right=lright;
# Line 822  DataLazy::~DataLazy() Line 821  DataLazy::~DataLazy()
821    For reasons of efficiency do not call this method on DataExpanded nodes.    For reasons of efficiency do not call this method on DataExpanded nodes.
822  */  */
823  DataReady_ptr  DataReady_ptr
824  DataLazy::collapseToReady()  DataLazy::collapseToReady() const
825  {  {
826    if (m_readytype=='E')    if (m_readytype=='E')
827    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
828      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
829    }    }
830    if (m_op==IDENTITY)    if (m_op==IDENTITY)
# Line 843  DataLazy::collapseToReady() Line 842  DataLazy::collapseToReady()
842    switch(m_op)    switch(m_op)
843    {    {
844      case ADD:      case ADD:
845      result=left+right;          result=left+right;
846      break;          break;
847      case SUB:            case SUB:          
848      result=left-right;          result=left-right;
849      break;          break;
850      case MUL:            case MUL:          
851      result=left*right;          result=left*right;
852      break;          break;
853      case DIV:            case DIV:          
854      result=left/right;          result=left/right;
855      break;          break;
856      case SIN:      case SIN:
857      result=left.sin();            result=left.sin();      
858      break;          break;
859      case COS:      case COS:
860      result=left.cos();          result=left.cos();
861      break;          break;
862      case TAN:      case TAN:
863      result=left.tan();          result=left.tan();
864      break;          break;
865      case ASIN:      case ASIN:
866      result=left.asin();          result=left.asin();
867      break;          break;
868      case ACOS:      case ACOS:
869      result=left.acos();          result=left.acos();
870      break;          break;
871      case ATAN:      case ATAN:
872      result=left.atan();          result=left.atan();
873      break;          break;
874      case SINH:      case SINH:
875      result=left.sinh();          result=left.sinh();
876      break;          break;
877      case COSH:      case COSH:
878      result=left.cosh();          result=left.cosh();
879      break;          break;
880      case TANH:      case TANH:
881      result=left.tanh();          result=left.tanh();
882      break;          break;
883      case ERF:      case ERF:
884      result=left.erf();          result=left.erf();
885      break;          break;
886     case ASINH:     case ASINH:
887      result=left.asinh();          result=left.asinh();
888      break;          break;
889     case ACOSH:     case ACOSH:
890      result=left.acosh();          result=left.acosh();
891      break;          break;
892     case ATANH:     case ATANH:
893      result=left.atanh();          result=left.atanh();
894      break;          break;
895      case LOG10:      case LOG10:
896      result=left.log10();          result=left.log10();
897      break;          break;
898      case LOG:      case LOG:
899      result=left.log();          result=left.log();
900      break;          break;
901      case SIGN:      case SIGN:
902      result=left.sign();          result=left.sign();
903      break;          break;
904      case ABS:      case ABS:
905      result=left.abs();          result=left.abs();
906      break;          break;
907      case NEG:      case NEG:
908      result=left.neg();          result=left.neg();
909      break;          break;
910      case POS:      case POS:
911      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
912      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
913      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
914      break;          break;
915      case EXP:      case EXP:
916      result=left.exp();          result=left.exp();
917      break;          break;
918      case SQRT:      case SQRT:
919      result=left.sqrt();          result=left.sqrt();
920      break;          break;
921      case RECIP:      case RECIP:
922      result=left.oneOver();          result=left.oneOver();
923      break;          break;
924      case GZ:      case GZ:
925      result=left.wherePositive();          result=left.wherePositive();
926      break;          break;
927      case LZ:      case LZ:
928      result=left.whereNegative();          result=left.whereNegative();
929      break;          break;
930      case GEZ:      case GEZ:
931      result=left.whereNonNegative();          result=left.whereNonNegative();
932      break;          break;
933      case LEZ:      case LEZ:
934      result=left.whereNonPositive();          result=left.whereNonPositive();
935      break;          break;
936      case NEZ:      case NEZ:
937      result=left.whereNonZero(m_tol);          result=left.whereNonZero(m_tol);
938      break;          break;
939      case EZ:      case EZ:
940      result=left.whereZero(m_tol);          result=left.whereZero(m_tol);
941      break;          break;
942      case SYM:      case SYM:
943      result=left.symmetric();          result=left.symmetric();
944      break;          break;
945      case NSYM:      case NSYM:
946      result=left.nonsymmetric();          result=left.nonsymmetric();
947      break;          break;
948      case PROD:      case PROD:
949      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);          result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
950      break;          break;
951      case TRANS:      case TRANS:
952      result=left.transpose(m_axis_offset);          result=left.transpose(m_axis_offset);
953      break;          break;
954      case TRACE:      case TRACE:
955      result=left.trace(m_axis_offset);          result=left.trace(m_axis_offset);
956      break;          break;
957      case SWAP:      case SWAP:
958      result=left.swapaxes(m_axis_offset, m_transpose);          result=left.swapaxes(m_axis_offset, m_transpose);
959      break;          break;
960      case MINVAL:      case MINVAL:
961      result=left.minval();          result=left.minval();
962      break;          break;
963      case MAXVAL:      case MAXVAL:
964      result=left.minval();          result=left.minval();
965      break;          break;
966      default:      default:
967      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)+".");
968    }    }
969    return result.borrowReadyPtr();    return result.borrowReadyPtr();
970  }  }
# Line 977  DataLazy::collapseToReady() Line 976  DataLazy::collapseToReady()
976     the purpose of using DataLazy in the first place).     the purpose of using DataLazy in the first place).
977  */  */
978  void  void
979  DataLazy::collapse()  DataLazy::collapse() const
980  {  {
981    if (m_op==IDENTITY)    if (m_op==IDENTITY)
982    {    {
983      return;          return;
984    }    }
985    if (m_readytype=='E')    if (m_readytype=='E')
986    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
987      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
988    }    }
989    m_id=collapseToReady();    m_id=collapseToReady();
# Line 997  DataLazy::collapse() Line 996  DataLazy::collapse()
996    
997    
998  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
999      for (int j=0;j<onumsteps;++j)\          for (int j=0;j<onumsteps;++j)\
1000      {\          {\
1001        for (int i=0;i<numsteps;++i,resultp+=resultStep) \            for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1002        { \            { \
1003  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1004  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1005           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \               tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1006  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1007           lroffset+=leftstep; \               lroffset+=leftstep; \
1008           rroffset+=rightstep; \               rroffset+=rightstep; \
1009        }\            }\
1010        lroffset+=oleftstep;\            lroffset+=oleftstep;\
1011        rroffset+=orightstep;\            rroffset+=orightstep;\
1012      }          }
1013    
1014    
1015  // The result will be stored in m_samples  // The result will be stored in m_samples
1016  // The return value is a pointer to the DataVector, offset is the offset within the return value  // The return value is a pointer to the DataVector, offset is the offset within the return value
1017  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1018  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
1019  {  {
1020  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1021      // 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
1022    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1023    {    {
1024      collapse();          collapse();
1025    }    }
1026    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
1027    {    {
1028      const ValueType& vec=m_id->getVectorRO();      const ValueType& vec=m_id->getVectorRO();
1029      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
# Line 1043  if (&x<stackend[omp_get_thread_num()]) Line 1042  if (&x<stackend[omp_get_thread_num()])
1042    }    }
1043    if (m_sampleids[tid]==sampleNo)    if (m_sampleids[tid]==sampleNo)
1044    {    {
1045      roffset=tid*m_samplesize;          roffset=tid*m_samplesize;
1046      return &(m_samples);        // sample is already resolved          return &(m_samples);            // sample is already resolved
1047    }    }
1048    m_sampleids[tid]=sampleNo;    m_sampleids[tid]=sampleNo;
1049    
# Line 1064  if (&x<stackend[omp_get_thread_num()]) Line 1063  if (&x<stackend[omp_get_thread_num()])
1063    }    }
1064  }  }
1065    
1066  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1067  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
1068  {  {
1069      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1070      // since we only have one argument we don't need to think about only          // since we only have one argument we don't need to think about only
1071      // processing single points.          // processing single points.
1072      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1073    if (m_readytype!='E')    if (m_readytype!='E')
1074    {    {
1075      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
# Line 1079  DataLazy::resolveNodeUnary(int tid, int Line 1078  DataLazy::resolveNodeUnary(int tid, int
1078    {    {
1079      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1080    }    }
1081    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1082    const double* left=&((*leftres)[roffset]);    const double* left=&((*leftres)[roffset]);
1083    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1084    double* result=&(m_samples[roffset]);    double* result=&(m_samples[roffset]);
1085    switch (m_op)    switch (m_op)
1086    {    {
1087      case SIN:        case SIN:  
1088      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
1089      break;          break;
1090      case COS:      case COS:
1091      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
1092      break;          break;
1093      case TAN:      case TAN:
1094      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
1095      break;          break;
1096      case ASIN:      case ASIN:
1097      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1098      break;          break;
1099      case ACOS:      case ACOS:
1100      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1101      break;          break;
1102      case ATAN:      case ATAN:
1103      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1104      break;          break;
1105      case SINH:      case SINH:
1106      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1107      break;          break;
1108      case COSH:      case COSH:
1109      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1110      break;          break;
1111      case TANH:      case TANH:
1112      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1113      break;          break;
1114      case ERF:      case ERF:
1115  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1116      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");          throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1117  #else  #else
1118      tensor_unary_operation(m_samplesize, left, result, ::erf);          tensor_unary_operation(m_samplesize, left, result, ::erf);
1119      break;          break;
1120  #endif  #endif
1121     case ASINH:     case ASINH:
1122  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1123      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);          tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1124  #else  #else
1125      tensor_unary_operation(m_samplesize, left, result, ::asinh);          tensor_unary_operation(m_samplesize, left, result, ::asinh);
1126  #endif    #endif  
1127      break;          break;
1128     case ACOSH:     case ACOSH:
1129  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1130      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);          tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1131  #else  #else
1132      tensor_unary_operation(m_samplesize, left, result, ::acosh);          tensor_unary_operation(m_samplesize, left, result, ::acosh);
1133  #endif    #endif  
1134      break;          break;
1135     case ATANH:     case ATANH:
1136  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1137      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);          tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1138  #else  #else
1139      tensor_unary_operation(m_samplesize, left, result, ::atanh);          tensor_unary_operation(m_samplesize, left, result, ::atanh);
1140  #endif    #endif  
1141      break;          break;
1142      case LOG10:      case LOG10:
1143      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1144      break;          break;
1145      case LOG:      case LOG:
1146      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1147      break;          break;
1148      case SIGN:      case SIGN:
1149      tensor_unary_operation(m_samplesize, left, result, escript::fsign);          tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1150      break;          break;
1151      case ABS:      case ABS:
1152      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1153      break;          break;
1154      case NEG:      case NEG:
1155      tensor_unary_operation(m_samplesize, left, result, negate<double>());          tensor_unary_operation(m_samplesize, left, result, negate<double>());
1156      break;          break;
1157      case POS:      case POS:
1158      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
1159      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
1160      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
1161      break;          break;
1162      case EXP:      case EXP:
1163      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1164      break;          break;
1165      case SQRT:      case SQRT:
1166      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1167      break;          break;
1168      case RECIP:      case RECIP:
1169      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));          tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1170      break;          break;
1171      case GZ:      case GZ:
1172      tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));          tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1173      break;          break;
1174      case LZ:      case LZ:
1175      tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));          tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1176      break;          break;
1177      case GEZ:      case GEZ:
1178      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));
1179      break;          break;
1180      case LEZ:      case LEZ:
1181      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));
1182      break;          break;
1183  // 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
1184      case NEZ:      case NEZ:
1185      tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));          tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1186      break;          break;
1187      case EZ:      case EZ:
1188      tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));          tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1189      break;          break;
1190    
1191      default:      default:
1192      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1193    }    }
1194    return &(m_samples);    return &(m_samples);
1195  }  }
1196    
1197    
1198  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1199  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1200  {  {
1201      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1202      // 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
1203      // processing single points.          // processing single points.
1204      // we will also know we won't get identity nodes          // we will also know we won't get identity nodes
1205    if (m_readytype!='E')    if (m_readytype!='E')
1206    {    {
1207      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
# Line 1212  DataLazy::resolveNodeReduction(int tid, Line 1211  DataLazy::resolveNodeReduction(int tid,
1211      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1212    }    }
1213    size_t loffset=0;    size_t loffset=0;
1214    const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);    const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1215    
1216    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1217    unsigned int ndpps=getNumDPPSample();    unsigned int ndpps=getNumDPPSample();
# Line 1221  DataLazy::resolveNodeReduction(int tid, Line 1220  DataLazy::resolveNodeReduction(int tid,
1220    switch (m_op)    switch (m_op)
1221    {    {
1222      case MINVAL:      case MINVAL:
1223      {          {
1224        for (unsigned int z=0;z<ndpps;++z)            for (unsigned int z=0;z<ndpps;++z)
1225        {            {
1226          FMin op;              FMin op;
1227          *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());
1228          loffset+=psize;              loffset+=psize;
1229          result++;              result++;
1230        }            }
1231      }          }
1232      break;          break;
1233      case MAXVAL:      case MAXVAL:
1234      {          {
1235        for (unsigned int z=0;z<ndpps;++z)            for (unsigned int z=0;z<ndpps;++z)
1236        {            {
1237        FMax op;            FMax op;
1238        *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);
1239        loffset+=psize;            loffset+=psize;
1240        result++;            result++;
1241        }            }
1242      }          }
1243      break;          break;
1244      default:      default:
1245      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1246    }    }
1247    return &(m_samples);    return &(m_samples);
1248  }  }
1249    
1250  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1251  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1252  {  {
1253      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1254      // 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
1255      // processing single points.          // processing single points.
1256    if (m_readytype!='E')    if (m_readytype!='E')
1257    {    {
1258      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
# Line 1272  DataLazy::resolveNodeNP1OUT(int tid, int Line 1271  DataLazy::resolveNodeNP1OUT(int tid, int
1271    switch (m_op)    switch (m_op)
1272    {    {
1273      case SYM:      case SYM:
1274      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1275      {          {
1276          DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1277          subroffset+=step;              subroffset+=step;
1278          offset+=step;              offset+=step;
1279      }          }
1280      break;          break;
1281      case NSYM:      case NSYM:
1282      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1283      {          {
1284          DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);              DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1285          subroffset+=step;              subroffset+=step;
1286          offset+=step;              offset+=step;
1287      }          }
1288      break;          break;
1289      default:      default:
1290      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1291    }    }
1292    return &m_samples;    return &m_samples;
1293  }  }
1294    
1295  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1296  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1297  {  {
1298      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1299      // 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
1300      // processing single points.          // processing single points.
1301    if (m_readytype!='E')    if (m_readytype!='E')
1302    {    {
1303      throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
# Line 1319  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1318  DataLazy::resolveNodeNP1OUT_P(int tid, i
1318    switch (m_op)    switch (m_op)
1319    {    {
1320      case TRACE:      case TRACE:
1321      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1322      {          {
1323              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);
1324          subroffset+=instep;              subroffset+=instep;
1325          offset+=outstep;              offset+=outstep;
1326      }          }
1327      break;          break;
1328      case TRANS:      case TRANS:
1329      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1330      {          {
1331              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);
1332          subroffset+=instep;              subroffset+=instep;
1333          offset+=outstep;              offset+=outstep;
1334      }          }
1335      break;          break;
1336      default:      default:
1337      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1338    }    }
1339    return &m_samples;    return &m_samples;
1340  }  }
1341    
1342    
1343  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1344  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1345  {  {
1346    if (m_readytype!='E')    if (m_readytype!='E')
1347    {    {
# Line 1364  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1363  DataLazy::resolveNodeNP1OUT_2P(int tid,
1363    switch (m_op)    switch (m_op)
1364    {    {
1365      case SWAP:      case SWAP:
1366      for (loop=0;loop<numsteps;++loop)          for (loop=0;loop<numsteps;++loop)
1367      {          {
1368              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);
1369          subroffset+=instep;              subroffset+=instep;
1370          offset+=outstep;              offset+=outstep;
1371      }          }
1372      break;          break;
1373      default:      default:
1374      throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1375    }    }
1376    return &m_samples;    return &m_samples;
1377  }  }
1378    
1379  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1380  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1381  {  {
1382    if (m_readytype!='E')    if (m_readytype!='E')
1383    {    {
# Line 1394  DataLazy::resolveNodeCondEval(int tid, i Line 1393  DataLazy::resolveNodeCondEval(int tid, i
1393    const ValueType* srcres=0;    const ValueType* srcres=0;
1394    if ((*maskres)[subroffset]>0)    if ((*maskres)[subroffset]>0)
1395    {    {
1396      srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1397    }    }
1398    else    else
1399    {    {
1400      srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1401    }    }
1402    
1403    // Now we need to copy the result    // Now we need to copy the result
# Line 1406  DataLazy::resolveNodeCondEval(int tid, i Line 1405  DataLazy::resolveNodeCondEval(int tid, i
1405    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1406    for (int i=0;i<m_samplesize;++i)    for (int i=0;i<m_samplesize;++i)
1407    {    {
1408      m_samples[roffset+i]=(*srcres)[subroffset+i];            m_samples[roffset+i]=(*srcres)[subroffset+i];  
1409    }    }
1410    
1411    return &m_samples;    return &m_samples;
# Line 1421  DataLazy::resolveNodeCondEval(int tid, i Line 1420  DataLazy::resolveNodeCondEval(int tid, i
1420  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1421  // For example, 2+Vector.  // For example, 2+Vector.
1422  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1423  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1424  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1425  {  {
1426  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1427    
1428    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
1429      // first work out which of the children are expanded          // first work out which of the children are expanded
1430    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1431    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1432    if (!leftExp && !rightExp)    if (!leftExp && !rightExp)
1433    {    {
1434      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'.");
1435    }    }
1436    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1437    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1438    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1439    {    {
1440      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.");
1441    }    }
1442    size_t leftsize=m_left->getNoValues();    size_t leftsize=m_left->getNoValues();
1443    size_t rightsize=m_right->getNoValues();    size_t rightsize=m_right->getNoValues();
1444    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
1445    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
1446    int rightstep=0;    int rightstep=0;
1447    int numsteps=0;       // total number of steps for the inner loop    int numsteps=0;               // total number of steps for the inner loop
1448    int oleftstep=0;  // the o variables refer to the outer loop    int oleftstep=0;      // the o variables refer to the outer loop
1449    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
1450    int onumsteps=1;    int onumsteps=1;
1451        
1452    bool LES=(leftExp && leftScalar); // Left is an expanded scalar    bool LES=(leftExp && leftScalar);     // Left is an expanded scalar
1453    bool RES=(rightExp && rightScalar);    bool RES=(rightExp && rightScalar);
1454    bool LS=(!leftExp && leftScalar); // left is a single scalar    bool LS=(!leftExp && leftScalar);     // left is a single scalar
1455    bool RS=(!rightExp && rightScalar);    bool RS=(!rightExp && rightScalar);
1456    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar    bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1457    bool RN=(!rightExp && !rightScalar);    bool RN=(!rightExp && !rightScalar);
1458    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar    bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1459    bool REN=(rightExp && !rightScalar);    bool REN=(rightExp && !rightScalar);
1460    
1461    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
1462    {    {
1463      chunksize=m_left->getNumDPPSample()*leftsize;          chunksize=m_left->getNumDPPSample()*leftsize;
1464      leftstep=0;          leftstep=0;
1465      rightstep=0;          rightstep=0;
1466      numsteps=1;          numsteps=1;
1467    }    }
1468    else if (LES || RES)    else if (LES || RES)
1469    {    {
1470      chunksize=1;          chunksize=1;
1471      if (LES)        // left is an expanded scalar          if (LES)                // left is an expanded scalar
1472      {          {
1473          if (RS)                  if (RS)
1474          {                  {
1475             leftstep=1;                     leftstep=1;
1476             rightstep=0;                     rightstep=0;
1477             numsteps=m_left->getNumDPPSample();                     numsteps=m_left->getNumDPPSample();
1478          }                  }
1479          else        // RN or REN                  else            // RN or REN
1480          {                  {
1481             leftstep=0;                     leftstep=0;
1482             oleftstep=1;                     oleftstep=1;
1483             rightstep=1;                     rightstep=1;
1484             orightstep=(RN ? -(int)rightsize : 0);                     orightstep=(RN ? -(int)rightsize : 0);
1485             numsteps=rightsize;                     numsteps=rightsize;
1486             onumsteps=m_left->getNumDPPSample();                     onumsteps=m_left->getNumDPPSample();
1487          }                  }
1488      }          }
1489      else        // right is an expanded scalar          else            // right is an expanded scalar
1490      {          {
1491          if (LS)                  if (LS)
1492          {                  {
1493             rightstep=1;                     rightstep=1;
1494             leftstep=0;                     leftstep=0;
1495             numsteps=m_right->getNumDPPSample();                     numsteps=m_right->getNumDPPSample();
1496          }                  }
1497          else                  else
1498          {                  {
1499             rightstep=0;                     rightstep=0;
1500             orightstep=1;                     orightstep=1;
1501             leftstep=1;                     leftstep=1;
1502             oleftstep=(LN ? -(int)leftsize : 0);                     oleftstep=(LN ? -(int)leftsize : 0);
1503             numsteps=leftsize;                     numsteps=leftsize;
1504             onumsteps=m_right->getNumDPPSample();                     onumsteps=m_right->getNumDPPSample();
1505          }                  }
1506      }          }
1507    }    }
1508    else  // this leaves (LEN, RS), (LEN, RN) and their transposes    else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1509    {    {
1510      if (LEN)    // and Right will be a single value          if (LEN)        // and Right will be a single value
1511      {          {
1512          chunksize=rightsize;                  chunksize=rightsize;
1513          leftstep=rightsize;                  leftstep=rightsize;
1514          rightstep=0;                  rightstep=0;
1515          numsteps=m_left->getNumDPPSample();                  numsteps=m_left->getNumDPPSample();
1516          if (RS)                  if (RS)
1517          {                  {
1518             numsteps*=leftsize;                     numsteps*=leftsize;
1519          }                  }
1520      }          }
1521      else    // REN          else    // REN
1522      {          {
1523          chunksize=leftsize;                  chunksize=leftsize;
1524          rightstep=leftsize;                  rightstep=leftsize;
1525          leftstep=0;                  leftstep=0;
1526          numsteps=m_right->getNumDPPSample();                  numsteps=m_right->getNumDPPSample();
1527          if (LS)                  if (LS)
1528          {                  {
1529             numsteps*=rightsize;                     numsteps*=rightsize;
1530          }                  }
1531      }          }
1532    }    }
1533    
1534    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1535      // Get the values of sub-expressions          // Get the values of sub-expressions
1536    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1537    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1538  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1539  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
# Line 1549  LAZYDEBUG(cout << "Right res["<< rroffse Line 1548  LAZYDEBUG(cout << "Right res["<< rroffse
1548    
1549    
1550    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1551    double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved    double* resultp=&(m_samples[roffset]);                // results are stored at the vector offset we received
1552    switch(m_op)    switch(m_op)
1553    {    {
1554      case ADD:      case ADD:
1555          PROC_OP(NO_ARG,plus<double>());          PROC_OP(NO_ARG,plus<double>());
1556      break;          break;
1557      case SUB:      case SUB:
1558      PROC_OP(NO_ARG,minus<double>());          PROC_OP(NO_ARG,minus<double>());
1559      break;          break;
1560      case MUL:      case MUL:
1561      PROC_OP(NO_ARG,multiplies<double>());          PROC_OP(NO_ARG,multiplies<double>());
1562      break;          break;
1563      case DIV:      case DIV:
1564      PROC_OP(NO_ARG,divides<double>());          PROC_OP(NO_ARG,divides<double>());
1565      break;          break;
1566      case POW:      case POW:
1567         PROC_OP(double (double,double),::pow);         PROC_OP(double (double,double),::pow);
1568      break;          break;
1569      default:      default:
1570      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1571    }    }
1572  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1573    return &m_samples;    return &m_samples;
# Line 1578  LAZYDEBUG(cout << "Result res[" << roffs Line 1577  LAZYDEBUG(cout << "Result res[" << roffs
1577  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1578  // have already been collapsed to IDENTITY. So we must have at least one expanded child.  // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1579  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1580  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1581  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1582  {  {
1583  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1584    
1585    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
1586      // first work out which of the children are expanded          // first work out which of the children are expanded
1587    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1588    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1589    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1590    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
1591    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1592    
1593    int resultStep=getNoValues();    int resultStep=getNoValues();
# Line 1611  LAZYDEBUG(cout << "m_samplesize=" << m_s Line 1610  LAZYDEBUG(cout << "m_samplesize=" << m_s
1610  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1611  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1612    
1613    double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(m_samples[offset]);         // results are stored at the vector offset we received
1614    switch(m_op)    switch(m_op)
1615    {    {
1616      case PROD:      case PROD:
1617      for (int i=0;i<steps;++i,resultp+=resultStep)          for (int i=0;i<steps;++i,resultp+=resultStep)
1618      {          {
1619            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1620            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1621    
1622  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1623  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1624    
1625            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);
1626    
1627        lroffset+=leftStep;            lroffset+=leftStep;
1628        rroffset+=rightStep;            rroffset+=rightStep;
1629      }          }
1630      break;          break;
1631      default:      default:
1632      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1633    }    }
1634    roffset=offset;    roffset=offset;
1635    return &m_samples;    return &m_samples;
1636  }  }
1637    
1638    
1639  const DataTypes::ValueType*  const DataTypes::RealVectorType*
1640  DataLazy::resolveSample(int sampleNo, size_t& roffset)  DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1641  {  {
1642  #ifdef _OPENMP  #ifdef _OPENMP
1643      int tid=omp_get_thread_num();          int tid=omp_get_thread_num();
1644  #else  #else
1645      int tid=0;          int tid=0;
1646  #endif  #endif
1647    
1648  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1649      stackstart[tid]=&tid;          stackstart[tid]=&tid;
1650      stackend[tid]=&tid;          stackend[tid]=&tid;
1651      const DataTypes::ValueType* r=resolveNodeSample(tid, sampleNo, roffset);          const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1652      size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];          size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1653      #pragma omp critical          #pragma omp critical
1654      if (d>maxstackuse)          if (d>maxstackuse)
1655      {          {
1656  cout << "Max resolve Stack use " << d << endl;  cout << "Max resolve Stack use " << d << endl;
1657          maxstackuse=d;                  maxstackuse=d;
1658      }          }
1659      return r;          return r;
1660  #else  #else
1661      return resolveNodeSample(tid, sampleNo, roffset);          return resolveNodeSample(tid, sampleNo, roffset);
1662  #endif  #endif
1663  }  }
1664    
# Line 1669  void Line 1668  void
1668  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
1669  {  {
1670     if (m_op==IDENTITY)     if (m_op==IDENTITY)
1671      return;          return;
1672     DataReady_ptr p=resolveNodeWorker();     DataReady_ptr p=resolveNodeWorker();
1673     makeIdentity(p);     makeIdentity(p);
1674  }  }
# Line 1706  DataLazy::resolveGroupWorker(std::vector Line 1705  DataLazy::resolveGroupWorker(std::vector
1705  {  {
1706    if (dats.empty())    if (dats.empty())
1707    {    {
1708      return;          return;
1709    }    }
1710    vector<DataLazy*> work;    vector<DataLazy*> work;
1711    FunctionSpace fs=dats[0]->getFunctionSpace();    FunctionSpace fs=dats[0]->getFunctionSpace();
1712    bool match=true;    bool match=true;
1713    for (int i=dats.size()-1;i>=0;--i)    for (int i=dats.size()-1;i>=0;--i)
1714    {    {
1715      if (dats[i]->m_readytype!='E')          if (dats[i]->m_readytype!='E')
1716      {          {
1717          dats[i]->collapse();                  dats[i]->collapse();
1718      }          }
1719      if (dats[i]->m_op!=IDENTITY)          if (dats[i]->m_op!=IDENTITY)
1720      {          {
1721          work.push_back(dats[i]);                  work.push_back(dats[i]);
1722          if (fs!=dats[i]->getFunctionSpace())                  if (fs!=dats[i]->getFunctionSpace())
1723          {                  {
1724              match=false;                          match=false;
1725          }                  }
1726      }          }
1727    }    }
1728    if (work.empty())    if (work.empty())
1729    {    {
1730      return;     // no work to do          return;         // no work to do
1731    }    }
1732    if (match)    // all functionspaces match.  Yes I realise this is overly strict    if (match)    // all functionspaces match.  Yes I realise this is overly strict
1733    {     // 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
1734          // all the other functionspaces match.                  // all the other functionspaces match.
1735      vector<DataExpanded*> dep;          vector<DataExpanded*> dep;
1736      vector<ValueType*> vecs;          vector<ValueType*> vecs;
1737      for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1738      {          {
1739          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())));
1740          vecs.push_back(&(dep[i]->getVectorRW()));                  vecs.push_back(&(dep[i]->getVectorRW()));
1741      }          }
1742      int totalsamples=work[0]->getNumSamples();          int totalsamples=work[0]->getNumSamples();
1743      const ValueType* res=0; // Storage for answer          const ValueType* res=0; // Storage for answer
1744      int sample;          int sample;
1745      #pragma omp parallel private(sample, res)          #pragma omp parallel private(sample, res)
1746      {          {
1747          size_t roffset=0;              size_t roffset=0;
1748          #pragma omp for schedule(static)              #pragma omp for schedule(static)
1749          for (sample=0;sample<totalsamples;++sample)              for (sample=0;sample<totalsamples;++sample)
1750          {              {
1751          roffset=0;                  roffset=0;
1752          int j;                  int j;
1753          for (j=work.size()-1;j>=0;--j)                  for (j=work.size()-1;j>=0;--j)
1754          {                  {
1755  #ifdef _OPENMP  #ifdef _OPENMP
1756                  res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);                      res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1757  #else  #else
1758                  res=work[j]->resolveNodeSample(0,sample,roffset);                      res=work[j]->resolveNodeSample(0,sample,roffset);
1759  #endif  #endif
1760                  DataVector::size_type outoffset=dep[j]->getPointOffset(sample,0);                      RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1761                  memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(DataVector::ElementType));                      memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
1762          }                  }
1763          }              }
1764      }          }
1765      // 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
1766      for (int i=work.size()-1;i>=0;--i)          for (int i=work.size()-1;i>=0;--i)
1767      {          {
1768          work[i]->makeIdentity(boost::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));              work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1769      }          }
1770    }    }
1771    else  // functionspaces do not match    else  // functionspaces do not match
1772    {    {
1773      for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1774      {          {
1775          work[i]->resolveToIdentity();                  work[i]->resolveToIdentity();
1776      }          }
1777    }    }
1778  }  }
1779    
# Line 1784  DataLazy::resolveGroupWorker(std::vector Line 1783  DataLazy::resolveGroupWorker(std::vector
1783  DataReady_ptr  DataReady_ptr
1784  DataLazy::resolveNodeWorker()  DataLazy::resolveNodeWorker()
1785  {  {
1786    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
1787    {    {
1788      collapse();      collapse();
1789    }    }
1790    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.
1791    {    {
1792      return m_id;      return m_id;
1793    }    }
1794      // 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'
1795    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1796    ValueType& resvec=result->getVectorRW();    ValueType& resvec=result->getVectorRW();
1797    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1798    
1799    int sample;    int sample;
1800    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1801    const ValueType* res=0;   // Storage for answer    const ValueType* res=0;       // Storage for answer
1802  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1803    #pragma omp parallel private(sample,res)    #pragma omp parallel private(sample,res)
1804    {    {
1805      size_t roffset=0;          size_t roffset=0;
1806  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1807      stackstart[omp_get_thread_num()]=&roffset;          stackstart[omp_get_thread_num()]=&roffset;
1808      stackend[omp_get_thread_num()]=&roffset;          stackend[omp_get_thread_num()]=&roffset;
1809  #endif  #endif
1810      #pragma omp for schedule(static)          #pragma omp for schedule(static)
1811      for (sample=0;sample<totalsamples;++sample)          for (sample=0;sample<totalsamples;++sample)
1812      {          {
1813          roffset=0;                  roffset=0;
1814  #ifdef _OPENMP  #ifdef _OPENMP
1815              res=resolveNodeSample(omp_get_thread_num(),sample,roffset);                  res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1816  #else  #else
1817              res=resolveNodeSample(0,sample,roffset);                  res=resolveNodeSample(0,sample,roffset);
1818  #endif  #endif
1819  LAZYDEBUG(cout << "Sample #" << sample << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1820  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1821              DataVector::size_type outoffset=result->getPointOffset(sample,0);                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1822              memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));                  memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1823      }          }
1824    }    }
1825  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1826    for (int i=0;i<getNumberOfThreads();++i)    for (int i=0;i<getNumberOfThreads();++i)
1827    {    {
1828      size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);          size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1829  //  cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;  //      cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1830      if (r>maxstackuse)          if (r>maxstackuse)
1831      {          {
1832          maxstackuse=r;                  maxstackuse=r;
1833      }          }
1834    }    }
1835    cout << "Max resolve Stack use=" << maxstackuse << endl;    cout << "Max resolve Stack use=" << maxstackuse << endl;
1836  #endif  #endif
# Line 1845  DataLazy::toString() const Line 1844  DataLazy::toString() const
1844    oss << "Lazy Data: [depth=" << m_height<< "] ";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1845    switch (escriptParams.getLAZY_STR_FMT())    switch (escriptParams.getLAZY_STR_FMT())
1846    {    {
1847    case 1:   // tree format    case 1:       // tree format
1848      oss << endl;          oss << endl;
1849      intoTreeString(oss,"");          intoTreeString(oss,"");
1850      break;          break;
1851    case 2:   // just the depth    case 2:       // just the depth
1852      break;          break;
1853    default:    default:
1854      intoString(oss);          intoString(oss);
1855      break;          break;
1856    }    }
1857    return oss.str();    return oss.str();
1858  }  }
# Line 1866  DataLazy::intoString(ostringstream& oss) Line 1865  DataLazy::intoString(ostringstream& oss)
1865    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1866    {    {
1867    case G_IDENTITY:    case G_IDENTITY:
1868      if (m_id->isExpanded())          if (m_id->isExpanded())
1869      {          {
1870         oss << "E";             oss << "E";
1871      }          }
1872      else if (m_id->isTagged())          else if (m_id->isTagged())
1873      {          {
1874        oss << "T";            oss << "T";
1875      }          }
1876      else if (m_id->isConstant())          else if (m_id->isConstant())
1877      {          {
1878        oss << "C";            oss << "C";
1879      }          }
1880      else          else
1881      {          {
1882        oss << "?";            oss << "?";
1883      }          }
1884      oss << '@' << m_id.get();          oss << '@' << m_id.get();
1885      break;          break;
1886    case G_BINARY:    case G_BINARY:
1887      oss << '(';          oss << '(';
1888      m_left->intoString(oss);          m_left->intoString(oss);
1889      oss << ' ' << opToString(m_op) << ' ';          oss << ' ' << opToString(m_op) << ' ';
1890      m_right->intoString(oss);          m_right->intoString(oss);
1891      oss << ')';          oss << ')';
1892      break;          break;
1893    case G_UNARY:    case G_UNARY:
1894    case G_UNARY_P:    case G_UNARY_P:
1895    case G_NP1OUT:    case G_NP1OUT:
1896    case G_NP1OUT_P:    case G_NP1OUT_P:
1897    case G_REDUCTION:    case G_REDUCTION:
1898      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1899      m_left->intoString(oss);          m_left->intoString(oss);
1900      oss << ')';          oss << ')';
1901      break;          break;
1902    case G_TENSORPROD:    case G_TENSORPROD:
1903      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1904      m_left->intoString(oss);          m_left->intoString(oss);
1905      oss << ", ";          oss << ", ";
1906      m_right->intoString(oss);          m_right->intoString(oss);
1907      oss << ')';          oss << ')';
1908      break;          break;
1909    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1910      oss << opToString(m_op) << '(';          oss << opToString(m_op) << '(';
1911      m_left->intoString(oss);          m_left->intoString(oss);
1912      oss << ", " << m_axis_offset << ", " << m_transpose;          oss << ", " << m_axis_offset << ", " << m_transpose;
1913      oss << ')';          oss << ')';
1914      break;          break;
1915    case G_CONDEVAL:    case G_CONDEVAL:
1916      oss << opToString(m_op)<< '(' ;          oss << opToString(m_op)<< '(' ;
1917      m_mask->intoString(oss);          m_mask->intoString(oss);
1918      oss << " ? ";          oss << " ? ";
1919      m_left->intoString(oss);          m_left->intoString(oss);
1920      oss << " : ";          oss << " : ";
1921      m_right->intoString(oss);          m_right->intoString(oss);
1922      oss << ')';          oss << ')';
1923      break;          break;
1924    default:    default:
1925      oss << "UNKNOWN";          oss << "UNKNOWN";
1926    }    }
1927  }  }
1928    
# Line 1935  DataLazy::intoTreeString(ostringstream& Line 1934  DataLazy::intoTreeString(ostringstream&
1934    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1935    {    {
1936    case G_IDENTITY:    case G_IDENTITY:
1937      if (m_id->isExpanded())          if (m_id->isExpanded())
1938      {          {
1939         oss << "E";             oss << "E";
1940      }          }
1941      else if (m_id->isTagged())          else if (m_id->isTagged())
1942      {          {
1943        oss << "T";            oss << "T";
1944      }          }
1945      else if (m_id->isConstant())          else if (m_id->isConstant())
1946      {          {
1947        oss << "C";            oss << "C";
1948      }          }
1949      else          else
1950      {          {
1951        oss << "?";            oss << "?";
1952      }          }
1953      oss << '@' << m_id.get() << endl;          oss << '@' << m_id.get() << endl;
1954      break;          break;
1955    case G_BINARY:    case G_BINARY:
1956      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1957      indent+='.';          indent+='.';
1958      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1959      m_right->intoTreeString(oss, indent);          m_right->intoTreeString(oss, indent);
1960      break;          break;
1961    case G_UNARY:    case G_UNARY:
1962    case G_UNARY_P:    case G_UNARY_P:
1963    case G_NP1OUT:    case G_NP1OUT:
1964    case G_NP1OUT_P:    case G_NP1OUT_P:
1965    case G_REDUCTION:    case G_REDUCTION:
1966      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1967      indent+='.';          indent+='.';
1968      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1969      break;          break;
1970    case G_TENSORPROD:    case G_TENSORPROD:
1971      oss << opToString(m_op) << endl;          oss << opToString(m_op) << endl;
1972      indent+='.';          indent+='.';
1973      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1974      m_right->intoTreeString(oss, indent);          m_right->intoTreeString(oss, indent);
1975      break;          break;
1976    case G_NP1OUT_2P:    case G_NP1OUT_2P:
1977      oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;          oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1978      indent+='.';          indent+='.';
1979      m_left->intoTreeString(oss, indent);          m_left->intoTreeString(oss, indent);
1980      break;          break;
1981    default:    default:
1982      oss << "UNKNOWN";          oss << "UNKNOWN";
1983    }    }
1984  }  }
1985    
1986    
1987  DataAbstract*  DataAbstract*
1988  DataLazy::deepCopy()  DataLazy::deepCopy() const
1989  {  {
1990    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1991    {    {
1992    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1993    case G_UNARY:    case G_UNARY:
1994    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1995    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);
1996    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);
1997    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);
1998    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);
1999    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);
2000    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);
2001    default:    default:
2002      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)+".");
2003    }    }
2004  }  }
2005    
# Line 2012  DataLazy::deepCopy() Line 2011  DataLazy::deepCopy()
2011  // or it could be some function of the lengths of the DataReady instances which  // or it could be some function of the lengths of the DataReady instances which
2012  // form part of the expression.  // form part of the expression.
2013  // Rather than have people making assumptions, I have disabled the method.  // Rather than have people making assumptions, I have disabled the method.
2014  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
2015  DataLazy::getLength() const  DataLazy::getLength() const
2016  {  {
2017    throw DataException("getLength() does not make sense for lazy data.");    throw DataException("getLength() does not make sense for lazy data.");
# Line 2027  DataLazy::getSlice(const DataTypes::Regi Line 2026  DataLazy::getSlice(const DataTypes::Regi
2026    
2027    
2028  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
2029  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
2030  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
2031                   int dataPointNo)                   int dataPointNo)
2032  {  {
2033    if (m_op==IDENTITY)    if (m_op==IDENTITY)
2034    {    {
2035      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
2036    }    }
2037    if (m_readytype!='E')    if (m_readytype!='E')
2038    {    {
2039      collapse();          collapse();
2040      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
2041    }    }
2042    // 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
2043    // so we only need to know which child to ask    // so we only need to know which child to ask
2044    if (m_left->m_readytype=='E')    if (m_left->m_readytype=='E')
2045    {    {
2046      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
2047    }    }
2048    else    else
2049    {    {
2050      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
2051    }    }
2052  }  }
2053    
2054  // To do this we need to rely on our child nodes  // To do this we need to rely on our child nodes
2055  DataTypes::ValueType::size_type  DataTypes::RealVectorType::size_type
2056  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
2057                   int dataPointNo) const                   int dataPointNo) const
2058  {  {
2059    if (m_op==IDENTITY)    if (m_op==IDENTITY)
2060    {    {
2061      return m_id->getPointOffset(sampleNo,dataPointNo);          return m_id->getPointOffset(sampleNo,dataPointNo);
2062    }    }
2063    if (m_readytype=='E')    if (m_readytype=='E')
2064    {    {
# Line 2067  DataLazy::getPointOffset(int sampleNo, Line 2066  DataLazy::getPointOffset(int sampleNo,
2066      // so we only need to know which child to ask      // so we only need to know which child to ask
2067      if (m_left->m_readytype=='E')      if (m_left->m_readytype=='E')
2068      {      {
2069      return m_left->getPointOffset(sampleNo,dataPointNo);          return m_left->getPointOffset(sampleNo,dataPointNo);
2070      }      }
2071      else      else
2072      {      {
2073      return m_right->getPointOffset(sampleNo,dataPointNo);          return m_right->getPointOffset(sampleNo,dataPointNo);
2074      }      }
2075    }    }
2076    if (m_readytype=='C')    if (m_readytype=='C')
2077    {    {
2078      return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter          return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
2079    }    }
2080    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).");
2081  }  }
# Line 2086  DataLazy::getPointOffset(int sampleNo, Line 2085  DataLazy::getPointOffset(int sampleNo,
2085  void  void
2086  DataLazy::setToZero()  DataLazy::setToZero()
2087  {  {
2088  //   DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::RealVectorType v(getNoValues(),0);
2089  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2090  //   m_op=IDENTITY;  //   m_op=IDENTITY;
2091  //   m_right.reset();    //   m_right.reset();  
# Line 2094  DataLazy::setToZero() Line 2093  DataLazy::setToZero()
2093  //   m_readytype='C';  //   m_readytype='C';
2094  //   m_buffsRequired=1;  //   m_buffsRequired=1;
2095    
2096    privdebug=privdebug;  // to stop the compiler complaining about unused privdebug    (void)privdebug;  // to stop the compiler complaining about unused privdebug
2097    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2098  }  }
2099    
2100  bool  bool
2101  DataLazy::actsExpanded() const  DataLazy::actsExpanded() const
2102  {  {
2103      return (m_readytype=='E');          return (m_readytype=='E');
2104  }  }
2105    
2106  }   // end namespace  } // end namespace
2107    

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

  ViewVC Help
Powered by ViewVC 1.1.26