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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26