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

Legend:
Removed from v.1888  
changed lines
  Added in v.2770

  ViewVC Help
Powered by ViewVC 1.1.26