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

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

  ViewVC Help
Powered by ViewVC 1.1.26