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

branches/schroedinger/escript/src/DataLazy.cpp revision 1910 by jfenwick, Thu Oct 23 03:05:28 2008 UTC trunk/escriptcore/src/DataLazy.cpp revision 6072 by jfenwick, Thu Mar 17 00:34:01 2016 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2016 by The University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * 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"
18    #include "Data.h"
19    #include "DataTypes.h"
20    #include "EscriptParams.h"
21    #include "FunctionSpace.h"
22    #include "UnaryFuncs.h"    // for escript::fsign
23    #include "Utils.h"
24    #include "DataMaths.h"
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
27  #include <netcdfcpp.h>  #include <netcdfcpp.h>
28  #endif  #endif
29  #ifdef PASO_MPI  
30  #include <mpi.h>  #include <iomanip> // for some fancy formatting in debug
31  #endif  
32  #ifdef _OPENMP  using namespace escript::DataTypes;
33  #include <omp.h>  
34  #endif  #define NO_ARG
35  #include "FunctionSpace.h"  
36  #include "DataTypes.h"  // #define LAZYDEBUG(X) if (privdebug){X;}
37  #include "Data.h"  #define LAZYDEBUG(X)
38  #include "UnaryFuncs.h"     // for escript::fsign  namespace
39    {
40    bool privdebug=false;
41    
42    #define ENABLEDEBUG privdebug=true;
43    #define DISABLEDEBUG privdebug=false;
44    }
45    
46    // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl;resolveToIdentity();}
47    
48    // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "SIZE LIMIT EXCEEDED " << m_height << endl;resolveToIdentity();}
49    
50    
51    #define SIZELIMIT if (m_height>escript::escriptParams.getTOO_MANY_LEVELS())  {if (escript::escriptParams.getLAZY_VERBOSE()){cerr << "SIZE LIMIT EXCEEDED height=" << m_height << endl;}resolveToIdentity();}
52    
53  /*  /*
54  How does DataLazy work?  How does DataLazy work?
# Line 38  A special operation, IDENTITY, stores an Line 61  A special operation, IDENTITY, stores an
61  This means that all "internal" nodes in the structure are instances of DataLazy.  This means that all "internal" nodes in the structure are instances of DataLazy.
62    
63  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...  Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
64  Note that IDENITY is not considered a unary operation.  Note that IDENTITY is not considered a unary operation.
65    
66  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).  I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
67  It must however form a DAG (directed acyclic graph).  It must however form a DAG (directed acyclic graph).
# Line 46  I will refer to individual DataLazy obje Line 69  I will refer to individual DataLazy obje
69    
70  Each node also stores:  Each node also stores:
71  - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was  - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
72      evaluated.          evaluated.
73  - m_length ~ how many values would be stored in the answer if the expression was evaluated.  - m_buffsrequired ~ the large number of samples which would need to be kept simultaneously in order to
74  - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to          evaluate the expression.
     evaluate the expression.  
75  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
76    
77  When a new node is created, the above values are computed based on the values in the child nodes.  When a new node is created, the above values are computed based on the values in the child nodes.
# Line 69  The convention that I use, is that the r Line 91  The convention that I use, is that the r
91    
92  For expressions which evaluate to Constant or Tagged, there is a different evaluation method.  For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
93  The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.  The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
94    
95    To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
96    1) Add to the ES_optype.
97    2) determine what opgroup your operation belongs to (X)
98    3) add a string for the op to the end of ES_opstrings
99    4) increase ES_opcount
100    5) add an entry (X) to opgroups
101    6) add an entry to the switch in collapseToReady
102    7) add an entry to resolveX
103  */  */
104    
105    
# Line 78  using namespace boost; Line 109  using namespace boost;
109  namespace escript  namespace escript
110  {  {
111    
 const std::string&  
 opToString(ES_optype op);  
   
112  namespace  namespace
113  {  {
114    
115    
116    // enabling this will print out when ever the maximum stacksize used by resolve increases
117    // it assumes _OPENMP is also in use
118    //#define LAZY_STACK_PROF
119    
120    
121    
122    #ifndef _OPENMP
123      #ifdef LAZY_STACK_PROF
124      #undef LAZY_STACK_PROF
125      #endif
126    #endif
127    
128    
129    #ifdef LAZY_STACK_PROF
130    std::vector<void*> stackstart(getNumberOfThreads());
131    std::vector<void*> stackend(getNumberOfThreads());
132    size_t maxstackuse=0;
133    #endif
134    
135  enum ES_opgroup  enum ES_opgroup
136  {  {
137     G_UNKNOWN,     G_UNKNOWN,
138     G_IDENTITY,     G_IDENTITY,
139     G_BINARY,        // pointwise operations with two arguments     G_BINARY,            // pointwise operations with two arguments
140     G_UNARY      // pointwise operations with one argument     G_UNARY,             // pointwise operations with one argument
141       G_UNARY_P,           // pointwise operations with one argument, requiring a parameter
142       G_NP1OUT,            // non-pointwise op with one output
143       G_NP1OUT_P,          // non-pointwise op with one output requiring a parameter
144       G_TENSORPROD,        // general tensor product
145       G_NP1OUT_2P,         // non-pointwise op with one output requiring two params
146       G_REDUCTION,         // non-pointwise unary op with a scalar output
147       G_CONDEVAL
148  };  };
149    
150    
151    
152    
153  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
154              "sin","cos","tan",                          "sin","cos","tan",
155              "asin","acos","atan","sinh","cosh","tanh","erf",                          "asin","acos","atan","sinh","cosh","tanh","erf",
156              "asinh","acosh","atanh",                          "asinh","acosh","atanh",
157              "log10","log","sign","abs","neg","pos","exp","sqrt",                          "log10","log","sign","abs","neg","pos","exp","sqrt",
158              "1/","where>0","where<0","where>=0","where<=0"};                          "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
159  int ES_opcount=32;                          "symmetric","nonsymmetric",
160                            "prod",
161                            "transpose", "trace",
162                            "swapaxes",
163                            "minval", "maxval",
164                            "condEval"};
165    int ES_opcount=44;
166  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
167              G_UNARY,G_UNARY,G_UNARY, //10                          G_UNARY,G_UNARY,G_UNARY, //10
168              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17                          G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 17
169              G_UNARY,G_UNARY,G_UNARY,                    // 20                          G_UNARY,G_UNARY,G_UNARY,                                        // 20
170              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28                          G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
171              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};                          G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,          // 35
172                            G_NP1OUT,G_NP1OUT,
173                            G_TENSORPROD,
174                            G_NP1OUT_P, G_NP1OUT_P,
175                            G_NP1OUT_2P,
176                            G_REDUCTION, G_REDUCTION,
177                            G_CONDEVAL};
178  inline  inline
179  ES_opgroup  ES_opgroup
180  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 119  getOpgroup(ES_optype op) Line 186  getOpgroup(ES_optype op)
186  FunctionSpace  FunctionSpace
187  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
188  {  {
189      // perhaps this should call interpolate and throw or something?          // perhaps this should call interpolate and throw or something?
190      // maybe we need an interpolate node -          // maybe we need an interpolate node -
191      // that way, if interpolate is required in any other op we can just throw a          // that way, if interpolate is required in any other op we can just throw a
192      // programming error exception.          // programming error exception.
193    
194      FunctionSpace l=left->getFunctionSpace();
195      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
196      {    if (l!=r)
197          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
198      }      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
199      return left->getFunctionSpace();      if (res==1)
200        {
201            return l;
202        }
203        if (res==-1)
204        {
205            return r;
206        }
207        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
208      }
209      return l;
210  }  }
211    
212  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
213    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
214  DataTypes::ShapeType  DataTypes::ShapeType
215  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
216  {  {
217      if (left->getShape()!=right->getShape())          if (left->getShape()!=right->getShape())
218      {          {
219        if (getOpgroup(op)!=G_BINARY)            if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
220        {            {
221          throw DataException("Shapes not the name - shapes must match for (point)binary operations.");                  throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
222        }            }
223        if (left->getRank()==0)   // we need to allow scalar * anything  
224        {            if (left->getRank()==0)       // we need to allow scalar * anything
225          return right->getShape();            {
226        }                  return right->getShape();
227        if (right->getRank()==0)            }
228        {            if (right->getRank()==0)
229          return left->getShape();            {
230        }                  return left->getShape();
231        throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");            }
232      }            throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
233      return left->getShape();          }
234            return left->getShape();
235  }  }
236    
237  // determine the number of points in the result of "left op right"  // return the shape for "op left"
238  size_t  
239  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataTypes::ShapeType
240    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
241  {  {
242     switch (getOpgroup(op))          switch(op)
243     {          {
244     case G_BINARY: return left->getLength();          case TRANS:
245     case G_UNARY: return left->getLength();             {                    // for the scoping of variables
246     default:                  const DataTypes::ShapeType& s=left->getShape();
247      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");                  DataTypes::ShapeType sh;
248     }                  int rank=left->getRank();
249                    if (axis_offset<0 || axis_offset>rank)
250                    {
251                stringstream e;
252                e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
253                throw DataException(e.str());
254            }
255            for (int i=0; i<rank; i++)
256                    {
257                       int index = (axis_offset+i)%rank;
258               sh.push_back(s[index]); // Append to new shape
259            }
260                    return sh;
261               }
262            break;
263            case TRACE:
264               {
265                    int rank=left->getRank();
266                    if (rank<2)
267                    {
268                       throw DataException("Trace can only be computed for objects with rank 2 or greater.");
269                    }
270                    if ((axis_offset>rank-2) || (axis_offset<0))
271                    {
272                       throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
273                    }
274                    if (rank==2)
275                    {
276                       return DataTypes::scalarShape;
277                    }
278                    else if (rank==3)
279                    {
280                       DataTypes::ShapeType sh;
281                       if (axis_offset==0)
282                       {
283                            sh.push_back(left->getShape()[2]);
284                       }
285                       else         // offset==1
286                       {
287                            sh.push_back(left->getShape()[0]);
288                       }
289                       return sh;
290                    }
291                    else if (rank==4)
292                    {
293                       DataTypes::ShapeType sh;
294                       const DataTypes::ShapeType& s=left->getShape();
295                       if (axis_offset==0)
296                       {
297                            sh.push_back(s[2]);
298                            sh.push_back(s[3]);
299                       }
300                       else if (axis_offset==1)
301                       {
302                            sh.push_back(s[0]);
303                            sh.push_back(s[3]);
304                       }
305                       else         // offset==2
306                       {
307                            sh.push_back(s[0]);
308                            sh.push_back(s[1]);
309                       }
310                       return sh;
311                    }
312                    else            // unknown rank
313                    {
314                       throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
315                    }
316               }
317            break;
318            default:
319            throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
320            }
321  }  }
322    
323  // determine the number of samples requires to evaluate an expression combining left and right  DataTypes::ShapeType
324  int  SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  
325  {  {
326     switch(getOpgroup(op))       // This code taken from the Data.cpp swapaxes() method
327     {       // Some of the checks are probably redundant here
328     case G_IDENTITY: return 1;       int axis0_tmp,axis1_tmp;
329     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);       const DataTypes::ShapeType& s=left->getShape();
330     case G_UNARY: return max(left->getBuffsRequired(),1);       DataTypes::ShapeType out_shape;
331     default:       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
332      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
333     }       int rank=left->getRank();
334         if (rank<2) {
335            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
336         }
337         if (axis0<0 || axis0>rank-1) {
338            stringstream e;
339            e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
340            throw DataException(e.str());
341         }
342         if (axis1<0 || axis1>rank-1) {
343            stringstream e;
344            e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
345            throw DataException(e.str());
346         }
347         if (axis0 == axis1) {
348             throw DataException("Error - Data::swapaxes: axis indices must be different.");
349         }
350         if (axis0 > axis1) {
351             axis0_tmp=axis1;
352             axis1_tmp=axis0;
353         } else {
354             axis0_tmp=axis0;
355             axis1_tmp=axis1;
356         }
357         for (int i=0; i<rank; i++) {
358           if (i == axis0_tmp) {
359              out_shape.push_back(s[axis1_tmp]);
360           } else if (i == axis1_tmp) {
361              out_shape.push_back(s[axis0_tmp]);
362           } else {
363              out_shape.push_back(s[i]);
364           }
365         }
366        return out_shape;
367  }  }
368    
369    
370  }   // end anonymous namespace  // determine the output shape for the general tensor product operation
371    // the additional parameters return information required later for the product
372    // the majority of this code is copy pasted from C_General_Tensor_Product
373    DataTypes::ShapeType
374    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
375    {
376            
377      // Get rank and shape of inputs
378      int rank0 = left->getRank();
379      int rank1 = right->getRank();
380      const DataTypes::ShapeType& shape0 = left->getShape();
381      const DataTypes::ShapeType& shape1 = right->getShape();
382    
383      // Prepare for the loops of the product and verify compatibility of shapes
384      int start0=0, start1=0;
385      if (transpose == 0)           {}
386      else if (transpose == 1)      { start0 = axis_offset; }
387      else if (transpose == 2)      { start1 = rank1-axis_offset; }
388      else                          { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
389    
390      if (rank0<axis_offset)
391      {
392            throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
393      }
394    
395      // Adjust the shapes for transpose
396      DataTypes::ShapeType tmpShape0(rank0);        // pre-sizing the vectors rather
397      DataTypes::ShapeType tmpShape1(rank1);        // than using push_back
398      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
399      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
400    
401      // Prepare for the loops of the product
402      SL=1, SM=1, SR=1;
403      for (int i=0; i<rank0-axis_offset; i++)       {
404        SL *= tmpShape0[i];
405      }
406      for (int i=rank0-axis_offset; i<rank0; i++)   {
407        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
408          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
409        }
410        SM *= tmpShape0[i];
411      }
412      for (int i=axis_offset; i<rank1; i++)         {
413        SR *= tmpShape1[i];
414      }
415    
416      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
417      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);      
418      {                     // block to limit the scope of out_index
419         int out_index=0;
420         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
421         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
422      }
423    
424      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
425      {
426         ostringstream os;
427         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
428         throw DataException(os.str());
429      }
430    
431      return shape2;
432    }
433    
434    }       // end anonymous namespace
435    
436    
437    
# Line 198  opToString(ES_optype op) Line 446  opToString(ES_optype op)
446    return ES_opstrings[op];    return ES_opstrings[op];
447  }  }
448    
449    void DataLazy::LazyNodeSetup()
450    {
451    #ifdef _OPENMP
452        int numthreads=omp_get_max_threads();
453        m_samples.resize(numthreads*m_samplesize);
454        m_sampleids=new int[numthreads];
455        for (int i=0;i<numthreads;++i)
456        {
457            m_sampleids[i]=-1;  
458        }
459    #else
460        m_samples.resize(m_samplesize);
461        m_sampleids=new int[1];
462        m_sampleids[0]=-1;
463    #endif  // _OPENMP
464    }
465    
466    
467    // Creates an identity node
468  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
469      : parent(p->getFunctionSpace(),p->getShape()),          : parent(p->getFunctionSpace(),p->getShape())
470      m_op(IDENTITY)          ,m_sampleids(0),
471            m_samples(1)
472  {  {
473     if (p->isLazy())     if (p->isLazy())
474     {     {
475      // I don't want identity of Lazy.          // I don't want identity of Lazy.
476      // Question: Why would that be so bad?          // Question: Why would that be so bad?
477      // Answer: We assume that the child of ID is something we can call getVector on          // Answer: We assume that the child of ID is something we can call getVector on
478      throw DataException("Programmer error - attempt to create identity from a DataLazy.");          throw DataException("Programmer error - attempt to create identity from a DataLazy.");
479     }     }
480     else     else
481     {     {
482      m_id=dynamic_pointer_cast<DataReady>(p);          p->makeLazyShared();
483      if(p->isConstant()) {m_readytype='C';}          DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
484      else if(p->isExpanded()) {m_readytype='E';}          makeIdentity(dr);
485      else if (p->isTagged()) {m_readytype='T';}  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}  
486     }     }
487     m_length=p->getLength();  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
 cout << "(1)Lazy created with " << m_samplesize << endl;  
488  }  }
489    
   
   
   
490  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
491      : parent(left->getFunctionSpace(),left->getShape()),          : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
492      m_op(op)          m_op(op),
493            m_axis_offset(0),
494            m_transpose(0),
495            m_SL(0), m_SM(0), m_SR(0)
496  {  {
497     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
498     {     {
499      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");          throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
500     }     }
501    
502     DataLazy_ptr lleft;     DataLazy_ptr lleft;
503     if (!left->isLazy())     if (!left->isLazy())
504     {     {
505      lleft=DataLazy_ptr(new DataLazy(left));          lleft=DataLazy_ptr(new DataLazy(left));
506     }     }
507     else     else
508     {     {
509      lleft=dynamic_pointer_cast<DataLazy>(left);          lleft=dynamic_pointer_cast<DataLazy>(left);
510     }     }
511     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
    m_length=left->getLength();  
512     m_left=lleft;     m_left=lleft;
    m_buffsRequired=1;  
513     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
514       m_children=m_left->m_children+1;
515       m_height=m_left->m_height+1;
516       LazyNodeSetup();
517       SIZELIMIT
518  }  }
519    
520    
521    // In this constructor we need to consider interpolation
522  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
523      : parent(resultFS(left,right,op), resultShape(left,right,op)),          : parent(resultFS(left,right,op), resultShape(left,right,op)),
524      m_op(op)          m_op(op),
525            m_SL(0), m_SM(0), m_SR(0)
526  {  {
527     if (getOpgroup(op)!=G_BINARY)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
528       if ((getOpgroup(op)!=G_BINARY))
529       {
530            throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
531       }
532    
533       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
534       {
535            FunctionSpace fs=getFunctionSpace();
536            Data ltemp(left);
537            Data tmp(ltemp,fs);
538            left=tmp.borrowDataPtr();
539       }
540       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
541     {     {
542      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");          Data tmp(Data(right),getFunctionSpace());
543            right=tmp.borrowDataPtr();
544    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
545     }     }
546     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     left->operandCheck(*right);
547    
548       if (left->isLazy())                  // the children need to be DataLazy. Wrap them in IDENTITY if required
549     {     {
550      m_left=dynamic_pointer_cast<DataLazy>(left);          m_left=dynamic_pointer_cast<DataLazy>(left);
551    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
552     }     }
553     else     else
554     {     {
555      m_left=DataLazy_ptr(new DataLazy(left));          m_left=DataLazy_ptr(new DataLazy(left));
556    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
557     }     }
558     if (right->isLazy())     if (right->isLazy())
559     {     {
560      m_right=dynamic_pointer_cast<DataLazy>(right);          m_right=dynamic_pointer_cast<DataLazy>(right);
561    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
562     }     }
563     else     else
564     {     {
565      m_right=DataLazy_ptr(new DataLazy(right));          m_right=DataLazy_ptr(new DataLazy(right));
566    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
567     }     }
568     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
569     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
570     if (lt=='E' || rt=='E')     if (lt=='E' || rt=='E')
571     {     {
572      m_readytype='E';          m_readytype='E';
573     }     }
574     else if (lt=='T' || rt=='T')     else if (lt=='T' || rt=='T')
575     {     {
576      m_readytype='T';          m_readytype='T';
577     }     }
578     else     else
579     {     {
580      m_readytype='C';          m_readytype='C';
581     }     }
582     m_length=resultLength(m_left,m_right,m_op);     m_samplesize=getNumDPPSample()*getNoValues();
583     m_samplesize=getNumDPPSample()*getNoValues();         m_children=m_left->m_children+m_right->m_children+2;
584     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
585  cout << "(3)Lazy created with " << m_samplesize << endl;     LazyNodeSetup();
586       SIZELIMIT
587    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
588  }  }
589    
590    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
591            : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
592            m_op(op),
593            m_axis_offset(axis_offset),
594            m_transpose(transpose)
595    {
596       if ((getOpgroup(op)!=G_TENSORPROD))
597       {
598            throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
599       }
600       if ((transpose>2) || (transpose<0))
601       {
602            throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
603       }
604       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
605       {
606            FunctionSpace fs=getFunctionSpace();
607            Data ltemp(left);
608            Data tmp(ltemp,fs);
609            left=tmp.borrowDataPtr();
610       }
611       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
612       {
613            Data tmp(Data(right),getFunctionSpace());
614            right=tmp.borrowDataPtr();
615       }
616    //    left->operandCheck(*right);
617    
618  DataLazy::~DataLazy()     if (left->isLazy())                  // the children need to be DataLazy. Wrap them in IDENTITY if required
619       {
620            m_left=dynamic_pointer_cast<DataLazy>(left);
621       }
622       else
623       {
624            m_left=DataLazy_ptr(new DataLazy(left));
625       }
626       if (right->isLazy())
627       {
628            m_right=dynamic_pointer_cast<DataLazy>(right);
629       }
630       else
631       {
632            m_right=DataLazy_ptr(new DataLazy(right));
633       }
634       char lt=m_left->m_readytype;
635       char rt=m_right->m_readytype;
636       if (lt=='E' || rt=='E')
637       {
638            m_readytype='E';
639       }
640       else if (lt=='T' || rt=='T')
641       {
642            m_readytype='T';
643       }
644       else
645       {
646            m_readytype='C';
647       }
648       m_samplesize=getNumDPPSample()*getNoValues();
649       m_children=m_left->m_children+m_right->m_children+2;
650       m_height=max(m_left->m_height,m_right->m_height)+1;
651       LazyNodeSetup();
652       SIZELIMIT
653    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
654    }
655    
656    
657    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
658            : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
659            m_op(op),
660            m_axis_offset(axis_offset),
661            m_transpose(0),
662            m_tol(0)
663  {  {
664       if ((getOpgroup(op)!=G_NP1OUT_P))
665       {
666            throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
667       }
668       DataLazy_ptr lleft;
669       if (!left->isLazy())
670       {
671            lleft=DataLazy_ptr(new DataLazy(left));
672       }
673       else
674       {
675            lleft=dynamic_pointer_cast<DataLazy>(left);
676       }
677       m_readytype=lleft->m_readytype;
678       m_left=lleft;
679       m_samplesize=getNumDPPSample()*getNoValues();
680       m_children=m_left->m_children+1;
681       m_height=m_left->m_height+1;
682       LazyNodeSetup();
683       SIZELIMIT
684    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
685    }
686    
687    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
688            : parent(left->getFunctionSpace(), left->getShape()),
689            m_op(op),
690            m_axis_offset(0),
691            m_transpose(0),
692            m_tol(tol)
693    {
694       if ((getOpgroup(op)!=G_UNARY_P))
695       {
696            throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
697       }
698       DataLazy_ptr lleft;
699       if (!left->isLazy())
700       {
701            lleft=DataLazy_ptr(new DataLazy(left));
702       }
703       else
704       {
705            lleft=dynamic_pointer_cast<DataLazy>(left);
706       }
707       m_readytype=lleft->m_readytype;
708       m_left=lleft;
709       m_samplesize=getNumDPPSample()*getNoValues();
710       m_children=m_left->m_children+1;
711       m_height=m_left->m_height+1;
712       LazyNodeSetup();
713       SIZELIMIT
714    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
715    }
716    
717    
718    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
719            : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
720            m_op(op),
721            m_axis_offset(axis0),
722            m_transpose(axis1),
723            m_tol(0)
724    {
725       if ((getOpgroup(op)!=G_NP1OUT_2P))
726       {
727            throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
728       }
729       DataLazy_ptr lleft;
730       if (!left->isLazy())
731       {
732            lleft=DataLazy_ptr(new DataLazy(left));
733       }
734       else
735       {
736            lleft=dynamic_pointer_cast<DataLazy>(left);
737       }
738       m_readytype=lleft->m_readytype;
739       m_left=lleft;
740       m_samplesize=getNumDPPSample()*getNoValues();
741       m_children=m_left->m_children+1;
742       m_height=m_left->m_height+1;
743       LazyNodeSetup();
744       SIZELIMIT
745    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
746    }
747    
748    
749    namespace
750    {
751    
752        inline int max3(int a, int b, int c)
753        {
754            int t=(a>b?a:b);
755            return (t>c?t:c);
756    
757        }
758    }
759    
760    DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
761            : parent(left->getFunctionSpace(), left->getShape()),
762            m_op(CONDEVAL),
763            m_axis_offset(0),
764            m_transpose(0),
765            m_tol(0)
766    {
767    
768       DataLazy_ptr lmask;
769       DataLazy_ptr lleft;
770       DataLazy_ptr lright;
771       if (!mask->isLazy())
772       {
773            lmask=DataLazy_ptr(new DataLazy(mask));
774       }
775       else
776       {
777            lmask=dynamic_pointer_cast<DataLazy>(mask);
778       }
779       if (!left->isLazy())
780       {
781            lleft=DataLazy_ptr(new DataLazy(left));
782       }
783       else
784       {
785            lleft=dynamic_pointer_cast<DataLazy>(left);
786       }
787       if (!right->isLazy())
788       {
789            lright=DataLazy_ptr(new DataLazy(right));
790       }
791       else
792       {
793            lright=dynamic_pointer_cast<DataLazy>(right);
794       }
795       m_readytype=lmask->m_readytype;
796       if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
797       {
798            throw DataException("Programmer Error - condEval arguments must have the same readytype");
799       }
800       m_left=lleft;
801       m_right=lright;
802       m_mask=lmask;
803       m_samplesize=getNumDPPSample()*getNoValues();
804       m_children=m_left->m_children+m_right->m_children+m_mask->m_children+1;
805       m_height=max3(m_left->m_height,m_right->m_height,m_mask->m_height)+1;
806       LazyNodeSetup();
807       SIZELIMIT
808    LAZYDEBUG(cout << "(8)Lazy created with " << m_samplesize << endl;)
809  }  }
810    
811    
812  int  
813  DataLazy::getBuffsRequired() const  DataLazy::~DataLazy()
814  {  {
815      return m_buffsRequired;     delete[] m_sampleids;
816  }  }
817    
818    
# Line 315  DataLazy::getBuffsRequired() const Line 822  DataLazy::getBuffsRequired() const
822    For reasons of efficiency do not call this method on DataExpanded nodes.    For reasons of efficiency do not call this method on DataExpanded nodes.
823  */  */
824  DataReady_ptr  DataReady_ptr
825  DataLazy::collapseToReady()  DataLazy::collapseToReady() const
826  {  {
827    if (m_readytype=='E')    if (m_readytype=='E')
828    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
829      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
830    }    }
831    if (m_op==IDENTITY)    if (m_op==IDENTITY)
# Line 328  DataLazy::collapseToReady() Line 835  DataLazy::collapseToReady()
835    DataReady_ptr pleft=m_left->collapseToReady();    DataReady_ptr pleft=m_left->collapseToReady();
836    Data left(pleft);    Data left(pleft);
837    Data right;    Data right;
838    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
839    {    {
840      right=Data(m_right->collapseToReady());      right=Data(m_right->collapseToReady());
841    }    }
# Line 336  DataLazy::collapseToReady() Line 843  DataLazy::collapseToReady()
843    switch(m_op)    switch(m_op)
844    {    {
845      case ADD:      case ADD:
846      result=left+right;          result=left+right;
847      break;          break;
848      case SUB:            case SUB:          
849      result=left-right;          result=left-right;
850      break;          break;
851      case MUL:            case MUL:          
852      result=left*right;          result=left*right;
853      break;          break;
854      case DIV:            case DIV:          
855      result=left/right;          result=left/right;
856      break;          break;
857      case SIN:      case SIN:
858      result=left.sin();            result=left.sin();      
859      break;          break;
860      case COS:      case COS:
861      result=left.cos();          result=left.cos();
862      break;          break;
863      case TAN:      case TAN:
864      result=left.tan();          result=left.tan();
865      break;          break;
866      case ASIN:      case ASIN:
867      result=left.asin();          result=left.asin();
868      break;          break;
869      case ACOS:      case ACOS:
870      result=left.acos();          result=left.acos();
871      break;          break;
872      case ATAN:      case ATAN:
873      result=left.atan();          result=left.atan();
874      break;          break;
875      case SINH:      case SINH:
876      result=left.sinh();          result=left.sinh();
877      break;          break;
878      case COSH:      case COSH:
879      result=left.cosh();          result=left.cosh();
880      break;          break;
881      case TANH:      case TANH:
882      result=left.tanh();          result=left.tanh();
883      break;          break;
884      case ERF:      case ERF:
885      result=left.erf();          result=left.erf();
886      break;          break;
887     case ASINH:     case ASINH:
888      result=left.asinh();          result=left.asinh();
889      break;          break;
890     case ACOSH:     case ACOSH:
891      result=left.acosh();          result=left.acosh();
892      break;          break;
893     case ATANH:     case ATANH:
894      result=left.atanh();          result=left.atanh();
895      break;          break;
896      case LOG10:      case LOG10:
897      result=left.log10();          result=left.log10();
898      break;          break;
899      case LOG:      case LOG:
900      result=left.log();          result=left.log();
901      break;          break;
902      case SIGN:      case SIGN:
903      result=left.sign();          result=left.sign();
904      break;          break;
905      case ABS:      case ABS:
906      result=left.abs();          result=left.abs();
907      break;          break;
908      case NEG:      case NEG:
909      result=left.neg();          result=left.neg();
910      break;          break;
911      case POS:      case POS:
912      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
913      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
914      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
915      break;          break;
916      case EXP:      case EXP:
917      result=left.exp();          result=left.exp();
918      break;          break;
919      case SQRT:      case SQRT:
920      result=left.sqrt();          result=left.sqrt();
921      break;          break;
922      case RECIP:      case RECIP:
923      result=left.oneOver();          result=left.oneOver();
924      break;          break;
925      case GZ:      case GZ:
926      result=left.wherePositive();          result=left.wherePositive();
927      break;          break;
928      case LZ:      case LZ:
929      result=left.whereNegative();          result=left.whereNegative();
930      break;          break;
931      case GEZ:      case GEZ:
932      result=left.whereNonNegative();          result=left.whereNonNegative();
933      break;          break;
934      case LEZ:      case LEZ:
935      result=left.whereNonPositive();          result=left.whereNonPositive();
936      break;          break;
937        case NEZ:
938            result=left.whereNonZero(m_tol);
939            break;
940        case EZ:
941            result=left.whereZero(m_tol);
942            break;
943        case SYM:
944            result=left.symmetric();
945            break;
946        case NSYM:
947            result=left.nonsymmetric();
948            break;
949        case PROD:
950            result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
951            break;
952        case TRANS:
953            result=left.transpose(m_axis_offset);
954            break;
955        case TRACE:
956            result=left.trace(m_axis_offset);
957            break;
958        case SWAP:
959            result=left.swapaxes(m_axis_offset, m_transpose);
960            break;
961        case MINVAL:
962            result=left.minval();
963            break;
964        case MAXVAL:
965            result=left.minval();
966            break;
967      default:      default:
968      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
969    }    }
970    return result.borrowReadyPtr();    return result.borrowReadyPtr();
971  }  }
# Line 440  DataLazy::collapseToReady() Line 977  DataLazy::collapseToReady()
977     the purpose of using DataLazy in the first place).     the purpose of using DataLazy in the first place).
978  */  */
979  void  void
980  DataLazy::collapse()  DataLazy::collapse() const
981  {  {
982    if (m_op==IDENTITY)    if (m_op==IDENTITY)
983    {    {
984      return;          return;
985    }    }
986    if (m_readytype=='E')    if (m_readytype=='E')
987    { // this is more an efficiency concern than anything else    {     // this is more an efficiency concern than anything else
988      throw DataException("Programmer Error - do not use collapse on Expanded data.");      throw DataException("Programmer Error - do not use collapse on Expanded data.");
989    }    }
990    m_id=collapseToReady();    m_id=collapseToReady();
991    m_op=IDENTITY;    m_op=IDENTITY;
992  }  }
993    
994  /*  // The result will be stored in m_samples
995    \brief Compute the value of the expression (binary operation) for the given sample.  // The return value is a pointer to the DataVector, offset is the offset within the return value
996    \return Vector which stores the value of the subexpression for the given sample.  const DataTypes::RealVectorType*
997    \param v A vector to store intermediate results.  DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
998    \param offset Index in v to begin storing results.  {
999    \param sampleNo Sample number to evaluate.  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1000    \param roffset (output parameter) the offset in the return vector where the result begins.          // collapse so we have a 'E' node or an IDENTITY for some other type
1001      if (m_readytype!='E' && m_op!=IDENTITY)
1002    The return value will be an existing vector so do not deallocate it.    {
1003    If the result is stored in v it should be stored at the offset given.          collapse();
1004    Everything from offset to the end of v should be considered available for this method to use.    }
1005  */    if (m_op==IDENTITY)  
1006  DataTypes::ValueType*    {
1007  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const      const RealVectorType& vec=m_id->getVectorRO();
1008        roffset=m_id->getPointOffset(sampleNo, 0);
1009    #ifdef LAZY_STACK_PROF
1010    int x;
1011    if (&x<stackend[omp_get_thread_num()])
1012    {
1013           stackend[omp_get_thread_num()]=&x;
1014    }
1015    #endif
1016        return &(vec);
1017      }
1018      if (m_readytype!='E')
1019      {
1020        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1021      }
1022      if (m_sampleids[tid]==sampleNo)
1023      {
1024            roffset=tid*m_samplesize;
1025            return &(m_samples);            // sample is already resolved
1026      }
1027      m_sampleids[tid]=sampleNo;
1028    
1029      switch (getOpgroup(m_op))
1030      {
1031      case G_UNARY:
1032      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1033      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1034      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1035      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1036      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1037      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1038      case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
1039      case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
1040      default:
1041        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1042      }
1043    }
1044    
1045    const DataTypes::RealVectorType*
1046    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
1047  {  {
1048      // we assume that any collapsing has been done before we get here          // we assume that any collapsing has been done before we get here
1049      // since we only have one argument we don't need to think about only          // since we only have one argument we don't need to think about only
1050      // processing single points.          // processing single points.
1051            // we will also know we won't get identity nodes
1052    if (m_readytype!='E')    if (m_readytype!='E')
1053    {    {
1054      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1055    }    }
1056    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);    if (m_op==IDENTITY)
1057    const double* left=&((*vleft)[roffset]);    {
1058    double* result=&(v[offset]);      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1059    roffset=offset;    }
1060      const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1061      const double* left=&((*leftres)[roffset]);
1062      roffset=m_samplesize*tid;
1063      double* result=&(m_samples[roffset]);
1064      escript::ESFunction operation=SINF;
1065    switch (m_op)    switch (m_op)
1066    {    {
1067      case SIN:        case SIN:
1068      tensor_unary_operation(m_samplesize, left, result, ::sin);      operation=SINF;
1069      break;      break;
1070      case COS:      case COS:
1071      tensor_unary_operation(m_samplesize, left, result, ::cos);          operation=COSF;
1072      break;      break;
1073      case TAN:      case TAN:
1074      tensor_unary_operation(m_samplesize, left, result, ::tan);          operation=TANF;
1075      break;      break;
1076      case ASIN:      case ASIN:
1077      tensor_unary_operation(m_samplesize, left, result, ::asin);          operation=ASINF;
1078      break;      break;
1079      case ACOS:      case ACOS:
1080      tensor_unary_operation(m_samplesize, left, result, ::acos);          operation=ACOSF;
1081      break;      break;
1082      case ATAN:      case ATAN:
1083      tensor_unary_operation(m_samplesize, left, result, ::atan);          operation=ATANF;
1084      break;      break;
1085      case SINH:      case SINH:
1086      tensor_unary_operation(m_samplesize, left, result, ::sinh);          operation=SINHF;
1087      break;      break;
1088      case COSH:      case COSH:
1089      tensor_unary_operation(m_samplesize, left, result, ::cosh);          operation=COSHF;
1090      break;      break;
1091      case TANH:      case TANH:
1092      tensor_unary_operation(m_samplesize, left, result, ::tanh);          operation=TANHF;
1093      break;      break;
1094      case ERF:      case ERF:
1095  #ifdef _WIN32          operation=ERFF;
     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::erf);  
1096      break;      break;
 #endif  
1097     case ASINH:     case ASINH:
1098  #ifdef _WIN32          operation=ASINHF;
     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 #endif    
1099      break;      break;
1100     case ACOSH:     case ACOSH:
1101  #ifdef _WIN32          operation=ACOSHF;
     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 #endif    
1102      break;      break;
1103     case ATANH:     case ATANH:
1104  #ifdef _WIN32          operation=ATANHF;
     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 #endif    
1105      break;      break;
1106      case LOG10:      case LOG10:
1107      tensor_unary_operation(m_samplesize, left, result, ::log10);          operation=LOG10F;
1108      break;      break;
1109      case LOG:      case LOG:
1110      tensor_unary_operation(m_samplesize, left, result, ::log);          operation=LOGF;
1111      break;      break;
1112      case SIGN:      case SIGN:
1113      tensor_unary_operation(m_samplesize, left, result, escript::fsign);          operation=SIGNF;
1114      break;      break;
1115      case ABS:      case ABS:
1116      tensor_unary_operation(m_samplesize, left, result, ::fabs);          operation=ABSF;
1117      break;      break;
1118      case NEG:      case NEG:
1119      tensor_unary_operation(m_samplesize, left, result, negate<double>());          operation=NEGF;
1120      break;      break;
1121      case POS:      case POS:
1122      // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
1123      // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
1124      throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
1125      break;          break;
1126      case EXP:      case EXP:
1127      tensor_unary_operation(m_samplesize, left, result, ::exp);          operation=EXPF;
1128      break;      break;
1129      case SQRT:      case SQRT:
1130      tensor_unary_operation(m_samplesize, left, result, ::sqrt);          operation=SQRTF;
1131      break;      break;
1132      case RECIP:      case RECIP:
1133      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));          operation=INVF;
1134      break;      break;
1135      case GZ:      case GZ:
1136      tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));          operation=GTZEROF;
1137      break;      break;
1138      case LZ:      case LZ:
1139      tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));          operation=LTZEROF;
1140      break;      break;
1141      case GEZ:      case GEZ:
1142      tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));          operation=GEZEROF;
1143      break;      break;
1144      case LEZ:      case LEZ:
1145      tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));          operation=LEZEROF;
1146        break;
1147    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1148        case NEZ:
1149            operation=NEQZEROF;
1150        break;
1151        case EZ:
1152            operation=EQZEROF;
1153      break;      break;
1154        default:
1155            throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1156      }
1157      tensor_unary_array_operation(m_samplesize,
1158                                 left,
1159                                 result,
1160                                 operation,
1161                                 m_tol);  
1162      return &(m_samples);
1163    }
1164    
1165    
1166    const DataTypes::RealVectorType*
1167    DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1168    {
1169            // we assume that any collapsing has been done before we get here
1170            // since we only have one argument we don't need to think about only
1171            // processing single points.
1172            // we will also know we won't get identity nodes
1173      if (m_readytype!='E')
1174      {
1175        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1176      }
1177      if (m_op==IDENTITY)
1178      {
1179        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1180      }
1181      size_t loffset=0;
1182      const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1183    
1184      roffset=m_samplesize*tid;
1185      unsigned int ndpps=getNumDPPSample();
1186      unsigned int psize=DataTypes::noValues(m_left->getShape());
1187      double* result=&(m_samples[roffset]);
1188      switch (m_op)
1189      {
1190        case MINVAL:
1191            {
1192              for (unsigned int z=0;z<ndpps;++z)
1193              {
1194                FMin op;
1195                *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1196                loffset+=psize;
1197                result++;
1198              }
1199            }
1200            break;
1201        case MAXVAL:
1202            {
1203              for (unsigned int z=0;z<ndpps;++z)
1204              {
1205              FMax op;
1206              *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1207              loffset+=psize;
1208              result++;
1209              }
1210            }
1211            break;
1212      default:      default:
1213      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1214    }    }
1215    return &v;    return &(m_samples);
1216  }  }
1217    
1218    const DataTypes::RealVectorType*
1219    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1220    {
1221            // we assume that any collapsing has been done before we get here
1222            // since we only have one argument we don't need to think about only
1223            // processing single points.
1224      if (m_readytype!='E')
1225      {
1226        throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1227      }
1228      if (m_op==IDENTITY)
1229      {
1230        throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1231      }
1232      size_t subroffset;
1233      const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1234      roffset=m_samplesize*tid;
1235      size_t loop=0;
1236      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1237      size_t step=getNoValues();
1238      size_t offset=roffset;
1239      switch (m_op)
1240      {
1241        case SYM:
1242            for (loop=0;loop<numsteps;++loop)
1243            {
1244                DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1245                subroffset+=step;
1246                offset+=step;
1247            }
1248            break;
1249        case NSYM:
1250            for (loop=0;loop<numsteps;++loop)
1251            {
1252                DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1253                subroffset+=step;
1254                offset+=step;
1255            }
1256            break;
1257        default:
1258            throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1259      }
1260      return &m_samples;
1261    }
1262    
1263    const DataTypes::RealVectorType*
1264    DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1265    {
1266            // we assume that any collapsing has been done before we get here
1267            // since we only have one argument we don't need to think about only
1268            // processing single points.
1269      if (m_readytype!='E')
1270      {
1271        throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1272      }
1273      if (m_op==IDENTITY)
1274      {
1275        throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1276      }
1277      size_t subroffset;
1278      size_t offset;
1279      const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1280      roffset=m_samplesize*tid;
1281      offset=roffset;
1282      size_t loop=0;
1283      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1284      size_t outstep=getNoValues();
1285      size_t instep=m_left->getNoValues();
1286      switch (m_op)
1287      {
1288        case TRACE:
1289            for (loop=0;loop<numsteps;++loop)
1290            {
1291                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1292                subroffset+=instep;
1293                offset+=outstep;
1294            }
1295            break;
1296        case TRANS:
1297            for (loop=0;loop<numsteps;++loop)
1298            {
1299                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1300                subroffset+=instep;
1301                offset+=outstep;
1302            }
1303            break;
1304        default:
1305            throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1306      }
1307      return &m_samples;
1308    }
1309    
1310    
1311    const DataTypes::RealVectorType*
1312    DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1313    {
1314      if (m_readytype!='E')
1315      {
1316        throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1317      }
1318      if (m_op==IDENTITY)
1319      {
1320        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1321      }
1322      size_t subroffset;
1323      size_t offset;
1324      const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1325      roffset=m_samplesize*tid;
1326      offset=roffset;
1327      size_t loop=0;
1328      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1329      size_t outstep=getNoValues();
1330      size_t instep=m_left->getNoValues();
1331      switch (m_op)
1332      {
1333        case SWAP:
1334            for (loop=0;loop<numsteps;++loop)
1335            {
1336                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1337                subroffset+=instep;
1338                offset+=outstep;
1339            }
1340            break;
1341        default:
1342            throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1343      }
1344      return &m_samples;
1345    }
1346    
1347  #define PROC_OP(X) \  const DataTypes::RealVectorType*
1348      for (int i=0;i<steps;++i,resultp+=resultStep) \  DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1349      { \  {
1350         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    if (m_readytype!='E')
1351         lroffset+=leftStep; \    {
1352         rroffset+=rightStep; \      throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1353      }    }
1354      if (m_op!=CONDEVAL)
1355      {
1356        throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1357      }
1358      size_t subroffset;
1359    
1360      const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1361      const RealVectorType* srcres=0;
1362      if ((*maskres)[subroffset]>0)
1363      {
1364            srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1365      }
1366      else
1367      {
1368            srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1369      }
1370    
1371      // Now we need to copy the result
1372    
1373      roffset=m_samplesize*tid;
1374      for (int i=0;i<m_samplesize;++i)
1375      {
1376            m_samples[roffset+i]=(*srcres)[subroffset+i];  
1377      }
1378    
1379      return &m_samples;
1380    }
1381    
 /*  
   \brief Compute the value of the expression (binary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
1382  // 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
1383  // 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.
1384  // If both children are expanded, then we can process them in a single operation (we treat  // If both children are expanded, then we can process them in a single operation (we treat
# Line 618  DataLazy::resolveUnary(ValueType& v, siz Line 1388  DataLazy::resolveUnary(ValueType& v, siz
1388  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1389  // For example, 2+Vector.  // For example, 2+Vector.
1390  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1391  DataTypes::ValueType*  const DataTypes::RealVectorType*
1392  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1393  {  {
1394  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1395    
1396    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
1397      // first work out which of the children are expanded          // first work out which of the children are expanded
1398    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1399    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1400    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1401    int steps=(bigloops?1:getNumDPPSample());    {
1402    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint          throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1403    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1404    {    bool leftScalar=(m_left->getRank()==0);
1405      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1406      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1407      chunksize=1;    // for scalar    {
1408    }              throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1409    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1410    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1411    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1412      // Get the values of sub-expressions    size_t chunksize=1;                   // how many doubles will be processed in one go
1413    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    int leftstep=0;               // how far should the left offset advance after each step
1414    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note    int rightstep=0;
1415      // the right child starts further along.    int numsteps=0;               // total number of steps for the inner loop
1416    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    int oleftstep=0;      // the o variables refer to the outer loop
1417      int orightstep=0;     // The outer loop is only required in cases where there is an extended scalar
1418      int onumsteps=1;
1419      
1420      bool LES=(leftExp && leftScalar);     // Left is an expanded scalar
1421      bool RES=(rightExp && rightScalar);
1422      bool LS=(!leftExp && leftScalar);     // left is a single scalar
1423      bool RS=(!rightExp && rightScalar);
1424      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1425      bool RN=(!rightExp && !rightScalar);
1426      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1427      bool REN=(rightExp && !rightScalar);
1428    
1429      if ((LES && RES) || (LEN && REN))     // both are Expanded scalars or both are expanded non-scalars
1430      {
1431            chunksize=m_left->getNumDPPSample()*leftsize;
1432            leftstep=0;
1433            rightstep=0;
1434            numsteps=1;
1435      }
1436      else if (LES || RES)
1437      {
1438            chunksize=1;
1439            if (LES)                // left is an expanded scalar
1440            {
1441                    if (RS)
1442                    {
1443                       leftstep=1;
1444                       rightstep=0;
1445                       numsteps=m_left->getNumDPPSample();
1446                    }
1447                    else            // RN or REN
1448                    {
1449                       leftstep=0;
1450                       oleftstep=1;
1451                       rightstep=1;
1452                       orightstep=(RN ? -(int)rightsize : 0);
1453                       numsteps=rightsize;
1454                       onumsteps=m_left->getNumDPPSample();
1455                    }
1456            }
1457            else            // right is an expanded scalar
1458            {
1459                    if (LS)
1460                    {
1461                       rightstep=1;
1462                       leftstep=0;
1463                       numsteps=m_right->getNumDPPSample();
1464                    }
1465                    else
1466                    {
1467                       rightstep=0;
1468                       orightstep=1;
1469                       leftstep=1;
1470                       oleftstep=(LN ? -(int)leftsize : 0);
1471                       numsteps=leftsize;
1472                       onumsteps=m_right->getNumDPPSample();
1473                    }
1474            }
1475      }
1476      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1477      {
1478            if (LEN)        // and Right will be a single value
1479            {
1480                    chunksize=rightsize;
1481                    leftstep=rightsize;
1482                    rightstep=0;
1483                    numsteps=m_left->getNumDPPSample();
1484                    if (RS)
1485                    {
1486                       numsteps*=leftsize;
1487                    }
1488            }
1489            else    // REN
1490            {
1491                    chunksize=leftsize;
1492                    rightstep=leftsize;
1493                    leftstep=0;
1494                    numsteps=m_right->getNumDPPSample();
1495                    if (LS)
1496                    {
1497                       numsteps*=rightsize;
1498                    }
1499            }
1500      }
1501    
1502      int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1503            // Get the values of sub-expressions
1504      const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1505      const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1506    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1507    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1508    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1509    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1510    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1511    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1512    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1513    
1514    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1515    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1516    
1517    
1518      roffset=m_samplesize*tid;
1519      double* resultp=&(m_samples[roffset]);                // results are stored at the vector offset we received
1520    switch(m_op)    switch(m_op)
1521    {    {
1522      case ADD:      case ADD:
1523      PROC_OP(plus<double>());          //PROC_OP(NO_ARG,plus<double>());
1524      break;        DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1525                 &(*left)[0],
1526                 &(*right)[0],
1527                 chunksize,
1528                 onumsteps,
1529                 numsteps,
1530                 resultStep,
1531                 leftstep,
1532                 rightstep,
1533                 oleftstep,
1534                 orightstep,
1535                 lroffset,
1536                 rroffset,
1537                 escript::ESFunction::PLUSF);  
1538            break;
1539      case SUB:      case SUB:
1540      PROC_OP(minus<double>());        DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1541      break;               &(*left)[0],
1542                 &(*right)[0],
1543                 chunksize,
1544                 onumsteps,
1545                 numsteps,
1546                 resultStep,
1547                 leftstep,
1548                 rightstep,
1549                 oleftstep,
1550                 orightstep,
1551                 lroffset,
1552                 rroffset,
1553                 escript::ESFunction::MINUSF);        
1554            //PROC_OP(NO_ARG,minus<double>());
1555            break;
1556      case MUL:      case MUL:
1557      PROC_OP(multiplies<double>());          //PROC_OP(NO_ARG,multiplies<double>());
1558      break;        DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1559                 &(*left)[0],
1560                 &(*right)[0],
1561                 chunksize,
1562                 onumsteps,
1563                 numsteps,
1564                 resultStep,
1565                 leftstep,
1566                 rightstep,
1567                 oleftstep,
1568                 orightstep,
1569                 lroffset,
1570                 rroffset,
1571                 escript::ESFunction::MULTIPLIESF);      
1572            break;
1573      case DIV:      case DIV:
1574      PROC_OP(divides<double>());          //PROC_OP(NO_ARG,divides<double>());
1575      break;        DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1576                 &(*left)[0],
1577                 &(*right)[0],
1578                 chunksize,
1579                 onumsteps,
1580                 numsteps,
1581                 resultStep,
1582                 leftstep,
1583                 rightstep,
1584                 oleftstep,
1585                 orightstep,
1586                 lroffset,
1587                 rroffset,
1588                 escript::ESFunction::DIVIDESF);          
1589            break;
1590      case POW:      case POW:
1591      PROC_OP(::pow);         //PROC_OP(double (double,double),::pow);
1592      break;        DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1593                 &(*left)[0],
1594                 &(*right)[0],
1595                 chunksize,
1596                 onumsteps,
1597                 numsteps,
1598                 resultStep,
1599                 leftstep,
1600                 rightstep,
1601                 oleftstep,
1602                 orightstep,
1603                 lroffset,
1604                 rroffset,
1605                 escript::ESFunction::POWF);          
1606            break;
1607      default:      default:
1608      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1609    }    }
1610    roffset=offset;    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1611    return &v;    return &m_samples;
1612  }  }
1613    
1614    
1615    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1616    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1617    // unlike the other resolve helpers, we must treat these datapoints separately.
1618    const DataTypes::RealVectorType*
1619    DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1620    {
1621    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1622    
1623  /*    size_t lroffset=0, rroffset=0;        // offsets in the left and right result vectors
1624    \brief Compute the value of the expression for the given sample.          // first work out which of the children are expanded
1625    \return Vector which stores the value of the subexpression for the given sample.    bool leftExp=(m_left->m_readytype=='E');
1626    \param v A vector to store intermediate results.    bool rightExp=(m_right->m_readytype=='E');
1627    \param offset Index in v to begin storing results.    int steps=getNumDPPSample();
1628    \param sampleNo Sample number to evaluate.    int leftStep=(leftExp? m_left->getNoValues() : 0);            // do not have scalars as input to this method
1629    \param roffset (output parameter) the offset in the return vector where the result begins.    int rightStep=(rightExp?m_right->getNoValues() : 0);
1630    
1631    The return value will be an existing vector so do not deallocate it.    int resultStep=getNoValues();
1632  */    roffset=m_samplesize*tid;
1633  // the vector and the offset are a place where the method could write its data if it wishes    size_t offset=roffset;
1634  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
1635  // Hence the return value to indicate where the data is actually stored.    const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1636  // Regardless, the storage should be assumed to be used, even if it isn't.  
1637      const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1638  // the roffset is the offset within the returned vector where the data begins  
1639  const DataTypes::ValueType*  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1640  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  cout << getNoValues() << endl;)
1641    
1642    
1643    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1644    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1645    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1646    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1647    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1648    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1649    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1650    
1651      double* resultp=&(m_samples[offset]);         // results are stored at the vector offset we received
1652      switch(m_op)
1653      {
1654        case PROD:
1655            for (int i=0;i<steps;++i,resultp+=resultStep)
1656            {
1657              const double *ptr_0 = &((*left)[lroffset]);
1658              const double *ptr_1 = &((*right)[rroffset]);
1659    
1660    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1661    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1662    
1663              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1664    
1665              lroffset+=leftStep;
1666              rroffset+=rightStep;
1667            }
1668            break;
1669        default:
1670            throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1671      }
1672      roffset=offset;
1673      return &m_samples;
1674    }
1675    
1676    
1677    const DataTypes::RealVectorType*
1678    DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1679  {  {
1680  cout << "Resolve sample " << toString() << endl;  #ifdef _OPENMP
1681      // collapse so we have a 'E' node or an IDENTITY for some other type          int tid=omp_get_thread_num();
1682    if (m_readytype!='E' && m_op!=IDENTITY)  #else
1683            int tid=0;
1684    #endif
1685    
1686    #ifdef LAZY_STACK_PROF
1687            stackstart[tid]=&tid;
1688            stackend[tid]=&tid;
1689            const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1690            size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1691            #pragma omp critical
1692            if (d>maxstackuse)
1693            {
1694    cout << "Max resolve Stack use " << d << endl;
1695                    maxstackuse=d;
1696            }
1697            return r;
1698    #else
1699            return resolveNodeSample(tid, sampleNo, roffset);
1700    #endif
1701    }
1702    
1703    
1704    // This needs to do the work of the identity constructor
1705    void
1706    DataLazy::resolveToIdentity()
1707    {
1708       if (m_op==IDENTITY)
1709            return;
1710       DataReady_ptr p=resolveNodeWorker();
1711       makeIdentity(p);
1712    }
1713    
1714    void DataLazy::makeIdentity(const DataReady_ptr& p)
1715    {
1716       m_op=IDENTITY;
1717       m_axis_offset=0;
1718       m_transpose=0;
1719       m_SL=m_SM=m_SR=0;
1720       m_children=m_height=0;
1721       m_id=p;
1722       if(p->isConstant()) {m_readytype='C';}
1723       else if(p->isExpanded()) {m_readytype='E';}
1724       else if (p->isTagged()) {m_readytype='T';}
1725       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1726       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1727       m_left.reset();
1728       m_right.reset();
1729    }
1730    
1731    
1732    DataReady_ptr
1733    DataLazy::resolve()
1734    {
1735        resolveToIdentity();
1736        return m_id;
1737    }
1738    
1739    
1740    /* This is really a static method but I think that caused problems in windows */
1741    void
1742    DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
1743    {
1744      if (dats.empty())
1745    {    {
1746      collapse();          return;
1747    }    }
1748    if (m_op==IDENTITY)      vector<DataLazy*> work;
1749      FunctionSpace fs=dats[0]->getFunctionSpace();
1750      bool match=true;
1751      for (int i=dats.size()-1;i>=0;--i)
1752    {    {
1753      const ValueType& vec=m_id->getVector();          if (dats[i]->m_readytype!='E')
1754      if (m_readytype=='C')          {
1755      {                  dats[i]->collapse();
1756      roffset=0;          }
1757      return &(vec);          if (dats[i]->m_op!=IDENTITY)
1758      }          {
1759      roffset=m_id->getPointOffset(sampleNo, 0);                  work.push_back(dats[i]);
1760      return &(vec);                  if (fs!=dats[i]->getFunctionSpace())
1761                    {
1762                            match=false;
1763                    }
1764            }
1765    }    }
1766    if (m_readytype!='E')    if (work.empty())
1767    {    {
1768      throw DataException("Programmer Error - Collapse did not produce an expanded node.");          return;         // no work to do
1769    }    }
1770    switch (getOpgroup(m_op))    if (match)    // all functionspaces match.  Yes I realise this is overly strict
1771      {             // it is possible that dats[0] is one of the objects which we discarded and
1772                    // all the other functionspaces match.
1773            vector<DataExpanded*> dep;
1774            vector<RealVectorType*> vecs;
1775            for (int i=0;i<work.size();++i)
1776            {
1777                    dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
1778                    vecs.push_back(&(dep[i]->getVectorRW()));
1779            }
1780            int totalsamples=work[0]->getNumSamples();
1781            const RealVectorType* res=0; // Storage for answer
1782            int sample;
1783            #pragma omp parallel private(sample, res)
1784            {
1785                size_t roffset=0;
1786                #pragma omp for schedule(static)
1787                for (sample=0;sample<totalsamples;++sample)
1788                {
1789                    roffset=0;
1790                    int j;
1791                    for (j=work.size()-1;j>=0;--j)
1792                    {
1793    #ifdef _OPENMP
1794                        res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1795    #else
1796                        res=work[j]->resolveNodeSample(0,sample,roffset);
1797    #endif
1798                        RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1799                        memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
1800                    }
1801                }
1802            }
1803            // Now we need to load the new results as identity ops into the lazy nodes
1804            for (int i=work.size()-1;i>=0;--i)
1805            {
1806                work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1807            }
1808      }
1809      else  // functionspaces do not match
1810    {    {
1811    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);          for (int i=0;i<work.size();++i)
1812    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);          {
1813    default:                  work[i]->resolveToIdentity();
1814      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");          }
1815    }    }
1816  }  }
1817    
1818    
1819  // To simplify the memory management, all threads operate on one large vector, rather than one each.  
1820  // Each sample is evaluated independently and copied into the result DataExpanded.  // This version of resolve uses storage in each node to hold results
1821  DataReady_ptr  DataReady_ptr
1822  DataLazy::resolve()  DataLazy::resolveNodeWorker()
1823  {  {
1824      if (m_readytype!='E')         // if the whole sub-expression is Constant or Tagged, then evaluate it normally
 cout << "Sample size=" << m_samplesize << endl;  
 cout << "Buffers=" << m_buffsRequired << endl;  
   
   if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally  
1825    {    {
1826      collapse();      collapse();
1827    }    }
1828    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.
1829    {    {
1830      return m_id;      return m_id;
1831    }    }
1832      // 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'
1833    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  RealVectorType(getNoValues()));
1834      // storage to evaluate its expression    RealVectorType& resvec=result->getVectorRW();
   int numthreads=1;  
 #ifdef _OPENMP  
   numthreads=getNumberOfThreads();  
   int threadnum=0;  
 #endif  
   ValueType v(numthreads*threadbuffersize);  
 cout << "Buffer created with size=" << v.size() << endl;  
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));  
   ValueType& resvec=result->getVector();  
1835    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1836    
1837    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1838    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1839    const ValueType* res=0;   // Vector storing the answer    const RealVectorType* res=0;       // Storage for answer
1840    size_t resoffset=0;       // where in the vector to find the answer  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1841    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)    #pragma omp parallel private(sample,res)
   for (sample=0;sample<totalsamples;++sample)  
1842    {    {
1843  cout << "################################# " << sample << endl;          size_t roffset=0;
1844    #ifdef LAZY_STACK_PROF
1845            stackstart[omp_get_thread_num()]=&roffset;
1846            stackend[omp_get_thread_num()]=&roffset;
1847    #endif
1848            #pragma omp for schedule(static)
1849            for (sample=0;sample<totalsamples;++sample)
1850            {
1851                    roffset=0;
1852  #ifdef _OPENMP  #ifdef _OPENMP
1853      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);                  res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1854  #else  #else
1855      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.                  res=resolveNodeSample(0,sample,roffset);
1856  #endif  #endif
1857  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1858      outoffset=result->getPointOffset(sample,0);  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1859  cerr << "offset=" << outoffset << endl;                  RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1860      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector                  memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1861      {          }
1862      resvec[outoffset]=(*res)[resoffset];    }
1863      }  #ifdef LAZY_STACK_PROF
1864  cerr << "*********************************" << endl;    for (int i=0;i<getNumberOfThreads();++i)
1865      {
1866            size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1867    //      cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1868            if (r>maxstackuse)
1869            {
1870                    maxstackuse=r;
1871            }
1872    }    }
1873      cout << "Max resolve Stack use=" << maxstackuse << endl;
1874    #endif
1875    return resptr;    return resptr;
1876  }  }
1877    
# Line 780  std::string Line 1879  std::string
1879  DataLazy::toString() const  DataLazy::toString() const
1880  {  {
1881    ostringstream oss;    ostringstream oss;
1882    oss << "Lazy Data:";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1883    intoString(oss);    switch (escriptParams.getLAZY_STR_FMT())
1884      {
1885      case 1:       // tree format
1886            oss << endl;
1887            intoTreeString(oss,"");
1888            break;
1889      case 2:       // just the depth
1890            break;
1891      default:
1892            intoString(oss);
1893            break;
1894      }
1895    return oss.str();    return oss.str();
1896  }  }
1897    
# Line 789  DataLazy::toString() const Line 1899  DataLazy::toString() const
1899  void  void
1900  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1901  {  {
1902    //    oss << "[" << m_children <<";"<<m_height <<"]";
1903    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1904    {    {
1905    case G_IDENTITY:    case G_IDENTITY:
1906      if (m_id->isExpanded())          if (m_id->isExpanded())
1907      {          {
1908         oss << "E";             oss << "E";
1909      }          }
1910      else if (m_id->isTagged())          else if (m_id->isTagged())
1911      {          {
1912        oss << "T";            oss << "T";
1913      }          }
1914      else if (m_id->isConstant())          else if (m_id->isConstant())
1915      {          {
1916        oss << "C";            oss << "C";
1917      }          }
1918      else          else
1919      {          {
1920        oss << "?";            oss << "?";
1921      }          }
1922      oss << '@' << m_id.get();          oss << '@' << m_id.get();
1923      break;          break;
1924    case G_BINARY:    case G_BINARY:
1925      oss << '(';          oss << '(';
1926      m_left->intoString(oss);          m_left->intoString(oss);
1927      oss << ' ' << opToString(m_op) << ' ';          oss << ' ' << opToString(m_op) << ' ';
1928      m_right->intoString(oss);          m_right->intoString(oss);
1929      oss << ')';          oss << ')';
1930      break;          break;
1931    case G_UNARY:    case G_UNARY:
1932      oss << opToString(m_op) << '(';    case G_UNARY_P:
1933      m_left->intoString(oss);    case G_NP1OUT:
1934      oss << ')';    case G_NP1OUT_P:
1935      break;    case G_REDUCTION:
1936            oss << opToString(m_op) << '(';
1937            m_left->intoString(oss);
1938            oss << ')';
1939            break;
1940      case G_TENSORPROD:
1941            oss << opToString(m_op) << '(';
1942            m_left->intoString(oss);
1943            oss << ", ";
1944            m_right->intoString(oss);
1945            oss << ')';
1946            break;
1947      case G_NP1OUT_2P:
1948            oss << opToString(m_op) << '(';
1949            m_left->intoString(oss);
1950            oss << ", " << m_axis_offset << ", " << m_transpose;
1951            oss << ')';
1952            break;
1953      case G_CONDEVAL:
1954            oss << opToString(m_op)<< '(' ;
1955            m_mask->intoString(oss);
1956            oss << " ? ";
1957            m_left->intoString(oss);
1958            oss << " : ";
1959            m_right->intoString(oss);
1960            oss << ')';
1961            break;
1962    default:    default:
1963      oss << "UNKNOWN";          oss << "UNKNOWN";
1964    }    }
1965  }  }
1966    
1967  // Note that in this case, deepCopy does not make copies of the leaves.  
1968  // Hopefully copy on write (or whatever we end up using) will take care of this.  void
1969    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1970    {
1971      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1972      switch (getOpgroup(m_op))
1973      {
1974      case G_IDENTITY:
1975            if (m_id->isExpanded())
1976            {
1977               oss << "E";
1978            }
1979            else if (m_id->isTagged())
1980            {
1981              oss << "T";
1982            }
1983            else if (m_id->isConstant())
1984            {
1985              oss << "C";
1986            }
1987            else
1988            {
1989              oss << "?";
1990            }
1991            oss << '@' << m_id.get() << endl;
1992            break;
1993      case G_BINARY:
1994            oss << opToString(m_op) << endl;
1995            indent+='.';
1996            m_left->intoTreeString(oss, indent);
1997            m_right->intoTreeString(oss, indent);
1998            break;
1999      case G_UNARY:
2000      case G_UNARY_P:
2001      case G_NP1OUT:
2002      case G_NP1OUT_P:
2003      case G_REDUCTION:
2004            oss << opToString(m_op) << endl;
2005            indent+='.';
2006            m_left->intoTreeString(oss, indent);
2007            break;
2008      case G_TENSORPROD:
2009            oss << opToString(m_op) << endl;
2010            indent+='.';
2011            m_left->intoTreeString(oss, indent);
2012            m_right->intoTreeString(oss, indent);
2013            break;
2014      case G_NP1OUT_2P:
2015            oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
2016            indent+='.';
2017            m_left->intoTreeString(oss, indent);
2018            break;
2019      default:
2020            oss << "UNKNOWN";
2021      }
2022    }
2023    
2024    
2025  DataAbstract*  DataAbstract*
2026  DataLazy::deepCopy()  DataLazy::deepCopy() const
2027  {  {
2028    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
2029    {    {
2030      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
2031      case G_UNARY:
2032      case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
2033      case G_UNARY_P:       return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
2034      case G_BINARY:        return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
2035      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
2036      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2037      case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
2038      case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2039      default:
2040            throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
2041    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
2042  }  }
2043    
2044    
2045  DataTypes::ValueType::size_type  
2046    // There is no single, natural interpretation of getLength on DataLazy.
2047    // Instances of DataReady can look at the size of their vectors.
2048    // For lazy though, it could be the size the data would be if it were resolved;
2049    // or it could be some function of the lengths of the DataReady instances which
2050    // form part of the expression.
2051    // Rather than have people making assumptions, I have disabled the method.
2052    DataTypes::RealVectorType::size_type
2053  DataLazy::getLength() const  DataLazy::getLength() const
2054  {  {
2055    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
2056  }  }
2057    
2058    
# Line 853  DataLazy::getSlice(const DataTypes::Regi Line 2062  DataLazy::getSlice(const DataTypes::Regi
2062    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
2063  }  }
2064    
2065  DataTypes::ValueType::size_type  
2066    // To do this we need to rely on our child nodes
2067    DataTypes::RealVectorType::size_type
2068    DataLazy::getPointOffset(int sampleNo,
2069                     int dataPointNo)
2070    {
2071      if (m_op==IDENTITY)
2072      {
2073            return m_id->getPointOffset(sampleNo,dataPointNo);
2074      }
2075      if (m_readytype!='E')
2076      {
2077            collapse();
2078            return m_id->getPointOffset(sampleNo,dataPointNo);
2079      }
2080      // at this point we do not have an identity node and the expression will be Expanded
2081      // so we only need to know which child to ask
2082      if (m_left->m_readytype=='E')
2083      {
2084            return m_left->getPointOffset(sampleNo,dataPointNo);
2085      }
2086      else
2087      {
2088            return m_right->getPointOffset(sampleNo,dataPointNo);
2089      }
2090    }
2091    
2092    // To do this we need to rely on our child nodes
2093    DataTypes::RealVectorType::size_type
2094  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
2095                   int dataPointNo) const                   int dataPointNo) const
2096  {  {
2097    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
2098      {
2099            return m_id->getPointOffset(sampleNo,dataPointNo);
2100      }
2101      if (m_readytype=='E')
2102      {
2103        // at this point we do not have an identity node and the expression will be Expanded
2104        // so we only need to know which child to ask
2105        if (m_left->m_readytype=='E')
2106        {
2107            return m_left->getPointOffset(sampleNo,dataPointNo);
2108        }
2109        else
2110        {
2111            return m_right->getPointOffset(sampleNo,dataPointNo);
2112        }
2113      }
2114      if (m_readytype=='C')
2115      {
2116            return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
2117      }
2118      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
2119  }  }
2120    
2121  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
2122  // to zero, all the tags from all the DataTags would be in the result.  // I have decided to let Data:: handle this issue.
 // However since they all have the same value (0) whether they are there or not should not matter.  
 // So I have decided that for all types this method will create a constant 0.  
 // It can be promoted up as required.  
 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management  
 // but we can deal with that if it arrises.  
2123  void  void
2124  DataLazy::setToZero()  DataLazy::setToZero()
2125  {  {
2126    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::RealVectorType v(getNoValues(),0);
2127    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2128    m_op=IDENTITY;  //   m_op=IDENTITY;
2129    m_right.reset();    //   m_right.reset();  
2130    m_left.reset();  //   m_left.reset();
2131    m_readytype='C';  //   m_readytype='C';
2132    m_buffsRequired=1;  //   m_buffsRequired=1;
2133    
2134      (void)privdebug;  // to stop the compiler complaining about unused privdebug
2135      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2136  }  }
2137    
2138  }   // end namespace  bool
2139    DataLazy::actsExpanded() const
2140    {
2141            return (m_readytype=='E');
2142    }
2143    
2144    } // end namespace
2145    

Legend:
Removed from v.1910  
changed lines
  Added in v.6072

  ViewVC Help
Powered by ViewVC 1.1.26