/[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_upto1946/escript/src/DataLazy.cpp revision 1950 by jfenwick, Thu Oct 30 00:59:34 2008 UTC trunk/escript/src/DataLazy.cpp revision 2271 by jfenwick, Mon Feb 16 05:08:29 2009 UTC
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31    #include <iomanip>      // for some fancy formatting in debug
32    
33    // #define LAZYDEBUG(X) if (privdebug){X;}
34    #define LAZYDEBUG(X)
35    namespace
36    {
37    bool privdebug=false;
38    
39    #define ENABLEDEBUG privdebug=true;
40    #define DISABLEDEBUG privdebug=false;
41    }
42    
43    // #define SIZELIMIT
44    // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
45    #define SIZELIMIT if ((m_height>7) || (m_children>127)) {resolveToIdentity();}
46    
47  /*  /*
48  How does DataLazy work?  How does DataLazy work?
49  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 48  I will refer to individual DataLazy obje Line 64  I will refer to individual DataLazy obje
64  Each node also stores:  Each node also stores:
65  - 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
66      evaluated.      evaluated.
 - m_length ~ how many values would be stored in the answer if the expression was evaluated.  
67  - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to  - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
68      evaluate the expression.      evaluate the expression.
69  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
# Line 70  The convention that I use, is that the r Line 85  The convention that I use, is that the r
85    
86  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.
87  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.
88    
89    To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
90    1) Add to the ES_optype.
91    2) determine what opgroup your operation belongs to (X)
92    3) add a string for the op to the end of ES_opstrings
93    4) increase ES_opcount
94    5) add an entry (X) to opgroups
95    6) add an entry to the switch in collapseToReady
96    7) add an entry to resolveX
97  */  */
98    
99    
# Line 79  using namespace boost; Line 103  using namespace boost;
103  namespace escript  namespace escript
104  {  {
105    
 const std::string&  
 opToString(ES_optype op);  
   
106  namespace  namespace
107  {  {
108    
# Line 90  enum ES_opgroup Line 111  enum ES_opgroup
111     G_UNKNOWN,     G_UNKNOWN,
112     G_IDENTITY,     G_IDENTITY,
113     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
114     G_UNARY      // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
115       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
116       G_NP1OUT,        // non-pointwise op with one output
117       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
118       G_TENSORPROD     // general tensor product
119  };  };
120    
121    
# Line 101  string ES_opstrings[]={"UNKNOWN","IDENTI Line 126  string ES_opstrings[]={"UNKNOWN","IDENTI
126              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
127              "asinh","acosh","atanh",              "asinh","acosh","atanh",
128              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
129              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
130  int ES_opcount=33;              "symmetric","nonsymmetric",
131                "prod",
132                "transpose", "trace"};
133    int ES_opcount=40;
134  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,
135              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
136              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
137              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
138              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
139              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
140                G_NP1OUT,G_NP1OUT,
141                G_TENSORPROD,
142                G_NP1OUT_P, G_NP1OUT_P};
143  inline  inline
144  ES_opgroup  ES_opgroup
145  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 143  resultFS(DataAbstract_ptr left, DataAbst Line 174  resultFS(DataAbstract_ptr left, DataAbst
174  }  }
175    
176  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
177    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
178  DataTypes::ShapeType  DataTypes::ShapeType
179  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
180  {  {
181      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
182      {      {
183        if (getOpgroup(op)!=G_BINARY)        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
184        {        {
185          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.");
186        }        }
# Line 165  resultShape(DataAbstract_ptr left, DataA Line 197  resultShape(DataAbstract_ptr left, DataA
197      return left->getShape();      return left->getShape();
198  }  }
199    
200  // determine the number of points in the result of "left op right"  // return the shape for "op left"
201  size_t  
202  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataTypes::ShapeType
203    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
204  {  {
205     switch (getOpgroup(op))      switch(op)
206     {      {
207     case G_BINARY: return left->getLength();          case TRANS:
208     case G_UNARY: return left->getLength();         {            // for the scoping of variables
209     default:          const DataTypes::ShapeType& s=left->getShape();
210      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");          DataTypes::ShapeType sh;
211     }          int rank=left->getRank();
212            if (axis_offset<0 || axis_offset>rank)
213            {
214                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
215                }
216                for (int i=0; i<rank; i++)
217            {
218               int index = (axis_offset+i)%rank;
219                   sh.push_back(s[index]); // Append to new shape
220                }
221            return sh;
222           }
223        break;
224        case TRACE:
225           {
226            int rank=left->getRank();
227            if (rank<2)
228            {
229               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
230            }
231            if ((axis_offset>rank-2) || (axis_offset<0))
232            {
233               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
234            }
235            if (rank==2)
236            {
237               return DataTypes::scalarShape;
238            }
239            else if (rank==3)
240            {
241               DataTypes::ShapeType sh;
242                   if (axis_offset==0)
243               {
244                    sh.push_back(left->getShape()[2]);
245                   }
246                   else     // offset==1
247               {
248                sh.push_back(left->getShape()[0]);
249                   }
250               return sh;
251            }
252            else if (rank==4)
253            {
254               DataTypes::ShapeType sh;
255               const DataTypes::ShapeType& s=left->getShape();
256                   if (axis_offset==0)
257               {
258                    sh.push_back(s[2]);
259                    sh.push_back(s[3]);
260                   }
261                   else if (axis_offset==1)
262               {
263                    sh.push_back(s[0]);
264                    sh.push_back(s[3]);
265                   }
266               else     // offset==2
267               {
268                sh.push_back(s[0]);
269                sh.push_back(s[1]);
270               }
271               return sh;
272            }
273            else        // unknown rank
274            {
275               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
276            }
277           }
278        break;
279            default:
280        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
281        }
282    }
283    
284    // determine the output shape for the general tensor product operation
285    // the additional parameters return information required later for the product
286    // the majority of this code is copy pasted from C_General_Tensor_Product
287    DataTypes::ShapeType
288    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
289    {
290        
291      // Get rank and shape of inputs
292      int rank0 = left->getRank();
293      int rank1 = right->getRank();
294      const DataTypes::ShapeType& shape0 = left->getShape();
295      const DataTypes::ShapeType& shape1 = right->getShape();
296    
297      // Prepare for the loops of the product and verify compatibility of shapes
298      int start0=0, start1=0;
299      if (transpose == 0)       {}
300      else if (transpose == 1)  { start0 = axis_offset; }
301      else if (transpose == 2)  { start1 = rank1-axis_offset; }
302      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
303    
304      if (rank0<axis_offset)
305      {
306        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
307      }
308    
309      // Adjust the shapes for transpose
310      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
311      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
312      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
313      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
314    
315      // Prepare for the loops of the product
316      SL=1, SM=1, SR=1;
317      for (int i=0; i<rank0-axis_offset; i++)   {
318        SL *= tmpShape0[i];
319      }
320      for (int i=rank0-axis_offset; i<rank0; i++)   {
321        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
322          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
323        }
324        SM *= tmpShape0[i];
325      }
326      for (int i=axis_offset; i<rank1; i++)     {
327        SR *= tmpShape1[i];
328      }
329    
330      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
331      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
332      {         // block to limit the scope of out_index
333         int out_index=0;
334         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
335         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
336      }
337    
338      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
339      {
340         ostringstream os;
341         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
342         throw DataException(os.str());
343      }
344    
345      return shape2;
346  }  }
347    
348  // determine the number of samples requires to evaluate an expression combining left and right  // determine the number of samples requires to evaluate an expression combining left and right
349    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
350    // The same goes for G_TENSORPROD
351    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
352    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
353    // multiple times
354  int  int
355  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
356  {  {
357     switch(getOpgroup(op))     switch(getOpgroup(op))
358     {     {
359     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
360     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
361     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY:
362       case G_UNARY_P: return max(left->getBuffsRequired(),1);
363       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
364       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
365       case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
366     default:     default:
367      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
368     }     }
# Line 210  opToString(ES_optype op) Line 386  opToString(ES_optype op)
386    
387    
388  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
389      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
     m_op(IDENTITY)  
390  {  {
391     if (p->isLazy())     if (p->isLazy())
392     {     {
# Line 222  DataLazy::DataLazy(DataAbstract_ptr p) Line 397  DataLazy::DataLazy(DataAbstract_ptr p)
397     }     }
398     else     else
399     {     {
400      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
401      if(p->isConstant()) {m_readytype='C';}      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
402      else if(p->isExpanded()) {m_readytype='E';}      makeIdentity(dr);
403      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.");}  
404     }     }
405     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;  
406  }  }
407    
408    
# Line 239  cout << "(1)Lazy created with " << m_sam Line 410  cout << "(1)Lazy created with " << m_sam
410    
411  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
412      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
413      m_op(op)      m_op(op),
414        m_axis_offset(0),
415        m_transpose(0),
416        m_SL(0), m_SM(0), m_SR(0)
417  {  {
418     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
419     {     {
420      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.");
421     }     }
422    
423     DataLazy_ptr lleft;     DataLazy_ptr lleft;
424     if (!left->isLazy())     if (!left->isLazy())
425     {     {
# Line 255  DataLazy::DataLazy(DataAbstract_ptr left Line 430  DataLazy::DataLazy(DataAbstract_ptr left
430      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
431     }     }
432     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
    m_length=left->getLength();  
433     m_left=lleft;     m_left=lleft;
434     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
435     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
436       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
437       m_children=m_left->m_children+1;
438       m_height=m_left->m_height+1;
439       SIZELIMIT
440  }  }
441    
442    
443  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
444  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
445      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
446      m_op(op)      m_op(op),
447        m_SL(0), m_SM(0), m_SR(0)
448  {  {
449     if (getOpgroup(op)!=G_BINARY)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
450       if ((getOpgroup(op)!=G_BINARY))
451     {     {
452      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.");
453     }     }
# Line 279  DataLazy::DataLazy(DataAbstract_ptr left Line 459  DataLazy::DataLazy(DataAbstract_ptr left
459      Data tmp(ltemp,fs);      Data tmp(ltemp,fs);
460      left=tmp.borrowDataPtr();      left=tmp.borrowDataPtr();
461     }     }
462     if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
463     {     {
464      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
465      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
466    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
467     }     }
468     left->operandCheck(*right);     left->operandCheck(*right);
469    
470     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
471     {     {
472      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
473    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
474       }
475       else
476       {
477        m_left=DataLazy_ptr(new DataLazy(left));
478    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
479       }
480       if (right->isLazy())
481       {
482        m_right=dynamic_pointer_cast<DataLazy>(right);
483    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
484       }
485       else
486       {
487        m_right=DataLazy_ptr(new DataLazy(right));
488    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
489       }
490       char lt=m_left->m_readytype;
491       char rt=m_right->m_readytype;
492       if (lt=='E' || rt=='E')
493       {
494        m_readytype='E';
495       }
496       else if (lt=='T' || rt=='T')
497       {
498        m_readytype='T';
499       }
500       else
501       {
502        m_readytype='C';
503       }
504       m_samplesize=getNumDPPSample()*getNoValues();
505       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
506       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
507       m_children=m_left->m_children+m_right->m_children+2;
508       m_height=max(m_left->m_height,m_right->m_height)+1;
509       SIZELIMIT
510    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
511    }
512    
513    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
514        : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
515        m_op(op),
516        m_axis_offset(axis_offset),
517        m_transpose(transpose)
518    {
519       if ((getOpgroup(op)!=G_TENSORPROD))
520       {
521        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
522       }
523       if ((transpose>2) || (transpose<0))
524       {
525        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
526       }
527       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
528       {
529        FunctionSpace fs=getFunctionSpace();
530        Data ltemp(left);
531        Data tmp(ltemp,fs);
532        left=tmp.borrowDataPtr();
533       }
534       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
535       {
536        Data tmp(Data(right),getFunctionSpace());
537        right=tmp.borrowDataPtr();
538       }
539    //    left->operandCheck(*right);
540    
541       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
542       {
543        m_left=dynamic_pointer_cast<DataLazy>(left);
544     }     }
545     else     else
546     {     {
# Line 316  DataLazy::DataLazy(DataAbstract_ptr left Line 568  DataLazy::DataLazy(DataAbstract_ptr left
568     {     {
569      m_readytype='C';      m_readytype='C';
570     }     }
571     m_length=resultLength(m_left,m_right,m_op);     m_samplesize=getNumDPPSample()*getNoValues();
572     m_samplesize=getNumDPPSample()*getNoValues();         m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
573     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
574  cout << "(3)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
575       m_height=max(m_left->m_height,m_right->m_height)+1;
576       SIZELIMIT
577    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
578    }
579    
580    
581    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
582        : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
583        m_op(op),
584        m_axis_offset(axis_offset),
585        m_transpose(0),
586        m_tol(0)
587    {
588       if ((getOpgroup(op)!=G_NP1OUT_P))
589       {
590        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
591       }
592       DataLazy_ptr lleft;
593       if (!left->isLazy())
594       {
595        lleft=DataLazy_ptr(new DataLazy(left));
596       }
597       else
598       {
599        lleft=dynamic_pointer_cast<DataLazy>(left);
600       }
601       m_readytype=lleft->m_readytype;
602       m_left=lleft;
603       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
604       m_samplesize=getNumDPPSample()*getNoValues();
605       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
606       m_children=m_left->m_children+1;
607       m_height=m_left->m_height+1;
608       SIZELIMIT
609    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
610  }  }
611    
612    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
613        : parent(left->getFunctionSpace(), left->getShape()),
614        m_op(op),
615        m_axis_offset(0),
616        m_transpose(0),
617        m_tol(tol)
618    {
619       if ((getOpgroup(op)!=G_UNARY_P))
620       {
621        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
622       }
623       DataLazy_ptr lleft;
624       if (!left->isLazy())
625       {
626        lleft=DataLazy_ptr(new DataLazy(left));
627       }
628       else
629       {
630        lleft=dynamic_pointer_cast<DataLazy>(left);
631       }
632       m_readytype=lleft->m_readytype;
633       m_left=lleft;
634       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
635       m_samplesize=getNumDPPSample()*getNoValues();
636       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
637       m_children=m_left->m_children+1;
638       m_height=m_left->m_height+1;
639       SIZELIMIT
640    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
641    }
642    
643  DataLazy::~DataLazy()  DataLazy::~DataLazy()
644  {  {
# Line 335  DataLazy::getBuffsRequired() const Line 652  DataLazy::getBuffsRequired() const
652  }  }
653    
654    
655    size_t
656    DataLazy::getMaxSampleSize() const
657    {
658        return m_maxsamplesize;
659    }
660    
661    
662    
663    size_t
664    DataLazy::getSampleBufferSize() const
665    {
666        return m_maxsamplesize*(max(1,m_buffsRequired));
667    }
668    
669  /*  /*
670    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
671    This does the work for the collapse method.    This does the work for the collapse method.
# Line 354  DataLazy::collapseToReady() Line 685  DataLazy::collapseToReady()
685    DataReady_ptr pleft=m_left->collapseToReady();    DataReady_ptr pleft=m_left->collapseToReady();
686    Data left(pleft);    Data left(pleft);
687    Data right;    Data right;
688    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
689    {    {
690      right=Data(m_right->collapseToReady());      right=Data(m_right->collapseToReady());
691    }    }
# Line 453  DataLazy::collapseToReady() Line 784  DataLazy::collapseToReady()
784      case LEZ:      case LEZ:
785      result=left.whereNonPositive();      result=left.whereNonPositive();
786      break;      break;
787        case NEZ:
788        result=left.whereNonZero(m_tol);
789        break;
790        case EZ:
791        result=left.whereZero(m_tol);
792        break;
793        case SYM:
794        result=left.symmetric();
795        break;
796        case NSYM:
797        result=left.nonsymmetric();
798        break;
799        case PROD:
800        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
801        break;
802        case TRANS:
803        result=left.transpose(m_axis_offset);
804        break;
805        case TRACE:
806        result=left.trace(m_axis_offset);
807        break;
808      default:      default:
809      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)+".");
810    }    }
811    return result.borrowReadyPtr();    return result.borrowReadyPtr();
812  }  }
# Line 481  DataLazy::collapse() Line 833  DataLazy::collapse()
833  }  }
834    
835  /*  /*
836    \brief Compute the value of the expression (binary operation) for the given sample.    \brief Compute the value of the expression (unary operation) for the given sample.
837    \return Vector which stores the value of the subexpression for the given sample.    \return Vector which stores the value of the subexpression for the given sample.
838    \param v A vector to store intermediate results.    \param v A vector to store intermediate results.
839    \param offset Index in v to begin storing results.    \param offset Index in v to begin storing results.
# Line 509  DataLazy::resolveUnary(ValueType& v, siz Line 861  DataLazy::resolveUnary(ValueType& v, siz
861    switch (m_op)    switch (m_op)
862    {    {
863      case SIN:        case SIN:  
864      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
865      break;      break;
866      case COS:      case COS:
867      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
868      break;      break;
869      case TAN:      case TAN:
870      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
871      break;      break;
872      case ASIN:      case ASIN:
873      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
874      break;      break;
875      case ACOS:      case ACOS:
876      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
877      break;      break;
878      case ATAN:      case ATAN:
879      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
880      break;      break;
881      case SINH:      case SINH:
882      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
883      break;      break;
884      case COSH:      case COSH:
885      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
886      break;      break;
887      case TANH:      case TANH:
888      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
889      break;      break;
890      case ERF:      case ERF:
891  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
892      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
893  #else  #else
894      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
895      break;      break;
896  #endif  #endif
897     case ASINH:     case ASINH:
898  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
899      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
900  #else  #else
901      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
902  #endif    #endif  
903      break;      break;
904     case ACOSH:     case ACOSH:
905  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
906      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
907  #else  #else
908      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
909  #endif    #endif  
910      break;      break;
911     case ATANH:     case ATANH:
912  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
913      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
914  #else  #else
915      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
916  #endif    #endif  
917      break;      break;
918      case LOG10:      case LOG10:
919      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
920      break;      break;
921      case LOG:      case LOG:
922      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
923      break;      break;
924      case SIGN:      case SIGN:
925      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
926      break;      break;
927      case ABS:      case ABS:
928      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
929      break;      break;
930      case NEG:      case NEG:
931      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 584  DataLazy::resolveUnary(ValueType& v, siz Line 936  DataLazy::resolveUnary(ValueType& v, siz
936      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
937      break;      break;
938      case EXP:      case EXP:
939      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
940      break;      break;
941      case SQRT:      case SQRT:
942      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
943      break;      break;
944      case RECIP:      case RECIP:
945      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 604  DataLazy::resolveUnary(ValueType& v, siz Line 956  DataLazy::resolveUnary(ValueType& v, siz
956      case LEZ:      case LEZ:
957      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));
958      break;      break;
959    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
960        case NEZ:
961        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
962        break;
963        case EZ:
964        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
965        break;
966    
967      default:      default:
968      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
# Line 615  DataLazy::resolveUnary(ValueType& v, siz Line 974  DataLazy::resolveUnary(ValueType& v, siz
974    
975    
976    
977  #define PROC_OP(X) \  
978      for (int i=0;i<steps;++i,resultp+=resultStep) \  /*
979      { \    \brief Compute the value of the expression (unary operation) for the given sample.
980         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    \return Vector which stores the value of the subexpression for the given sample.
981         lroffset+=leftStep; \    \param v A vector to store intermediate results.
982         rroffset+=rightStep; \    \param offset Index in v to begin storing results.
983      \param sampleNo Sample number to evaluate.
984      \param roffset (output parameter) the offset in the return vector where the result begins.
985    
986      The return value will be an existing vector so do not deallocate it.
987      If the result is stored in v it should be stored at the offset given.
988      Everything from offset to the end of v should be considered available for this method to use.
989    */
990    DataTypes::ValueType*
991    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
992    {
993        // we assume that any collapsing has been done before we get here
994        // since we only have one argument we don't need to think about only
995        // processing single points.
996      if (m_readytype!='E')
997      {
998        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
999      }
1000        // since we can't write the result over the input, we need a result offset further along
1001      size_t subroffset=roffset+m_samplesize;
1002    LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
1003      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1004      roffset=offset;
1005      size_t loop=0;
1006      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1007      size_t step=getNoValues();
1008      switch (m_op)
1009      {
1010        case SYM:
1011        for (loop=0;loop<numsteps;++loop)
1012        {
1013            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1014            subroffset+=step;
1015            offset+=step;
1016        }
1017        break;
1018        case NSYM:
1019        for (loop=0;loop<numsteps;++loop)
1020        {
1021            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1022            subroffset+=step;
1023            offset+=step;
1024        }
1025        break;
1026        default:
1027        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1028      }
1029      return &v;
1030    }
1031    
1032    /*
1033      \brief Compute the value of the expression (unary operation) for the given sample.
1034      \return Vector which stores the value of the subexpression for the given sample.
1035      \param v A vector to store intermediate results.
1036      \param offset Index in v to begin storing results.
1037      \param sampleNo Sample number to evaluate.
1038      \param roffset (output parameter) the offset in the return vector where the result begins.
1039    
1040      The return value will be an existing vector so do not deallocate it.
1041      If the result is stored in v it should be stored at the offset given.
1042      Everything from offset to the end of v should be considered available for this method to use.
1043    */
1044    DataTypes::ValueType*
1045    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1046    {
1047        // we assume that any collapsing has been done before we get here
1048        // since we only have one argument we don't need to think about only
1049        // processing single points.
1050      if (m_readytype!='E')
1051      {
1052        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1053      }
1054        // since we can't write the result over the input, we need a result offset further along
1055      size_t subroffset;
1056      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1057    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1058    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1059      roffset=offset;
1060      size_t loop=0;
1061      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1062      size_t outstep=getNoValues();
1063      size_t instep=m_left->getNoValues();
1064    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1065      switch (m_op)
1066      {
1067        case TRACE:
1068        for (loop=0;loop<numsteps;++loop)
1069        {
1070    size_t zz=sampleNo*getNumDPPSample()+loop;
1071    if (zz==5800)
1072    {
1073    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
1074    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1075    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1076    LAZYDEBUG(cerr << subroffset << endl;)
1077    LAZYDEBUG(cerr << "output=" << offset << endl;)
1078    }
1079                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1080    if (zz==5800)
1081    {
1082    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1083    }
1084            subroffset+=instep;
1085            offset+=outstep;
1086        }
1087        break;
1088        case TRANS:
1089        for (loop=0;loop<numsteps;++loop)
1090        {
1091                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1092            subroffset+=instep;
1093            offset+=outstep;
1094        }
1095        break;
1096        default:
1097        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1098      }
1099      return &v;
1100    }
1101    
1102    
1103    #define PROC_OP(TYPE,X)                               \
1104        for (int j=0;j<onumsteps;++j)\
1105        {\
1106          for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1107          { \
1108    LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1109    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1110             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1111    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1112             lroffset+=leftstep; \
1113             rroffset+=rightstep; \
1114          }\
1115          lroffset+=oleftstep;\
1116          rroffset+=orightstep;\
1117      }      }
1118    
1119  /*  /*
# Line 647  DataLazy::resolveUnary(ValueType& v, siz Line 1140  DataLazy::resolveUnary(ValueType& v, siz
1140  DataTypes::ValueType*  DataTypes::ValueType*
1141  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1142  {  {
1143  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1144    
1145    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
1146      // first work out which of the children are expanded      // first work out which of the children are expanded
1147    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1148    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1149    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1150    int steps=(bigloops?1:getNumDPPSample());    {
1151    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'.");
1152    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1153    {    bool leftScalar=(m_left->getRank()==0);
1154      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1155      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1156      chunksize=1;    // for scalar    {
1157    }          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1158    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1159    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1160    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1161      size_t chunksize=1;           // how many doubles will be processed in one go
1162      int leftstep=0;       // how far should the left offset advance after each step
1163      int rightstep=0;
1164      int numsteps=0;       // total number of steps for the inner loop
1165      int oleftstep=0;  // the o variables refer to the outer loop
1166      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1167      int onumsteps=1;
1168      
1169      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1170      bool RES=(rightExp && rightScalar);
1171      bool LS=(!leftExp && leftScalar); // left is a single scalar
1172      bool RS=(!rightExp && rightScalar);
1173      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1174      bool RN=(!rightExp && !rightScalar);
1175      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1176      bool REN=(rightExp && !rightScalar);
1177    
1178      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1179      {
1180        chunksize=m_left->getNumDPPSample()*leftsize;
1181        leftstep=0;
1182        rightstep=0;
1183        numsteps=1;
1184      }
1185      else if (LES || RES)
1186      {
1187        chunksize=1;
1188        if (LES)        // left is an expanded scalar
1189        {
1190            if (RS)
1191            {
1192               leftstep=1;
1193               rightstep=0;
1194               numsteps=m_left->getNumDPPSample();
1195            }
1196            else        // RN or REN
1197            {
1198               leftstep=0;
1199               oleftstep=1;
1200               rightstep=1;
1201               orightstep=(RN ? -(int)rightsize : 0);
1202               numsteps=rightsize;
1203               onumsteps=m_left->getNumDPPSample();
1204            }
1205        }
1206        else        // right is an expanded scalar
1207        {
1208            if (LS)
1209            {
1210               rightstep=1;
1211               leftstep=0;
1212               numsteps=m_right->getNumDPPSample();
1213            }
1214            else
1215            {
1216               rightstep=0;
1217               orightstep=1;
1218               leftstep=1;
1219               oleftstep=(LN ? -(int)leftsize : 0);
1220               numsteps=leftsize;
1221               onumsteps=m_right->getNumDPPSample();
1222            }
1223        }
1224      }
1225      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1226      {
1227        if (LEN)    // and Right will be a single value
1228        {
1229            chunksize=rightsize;
1230            leftstep=rightsize;
1231            rightstep=0;
1232            numsteps=m_left->getNumDPPSample();
1233            if (RS)
1234            {
1235               numsteps*=leftsize;
1236            }
1237        }
1238        else    // REN
1239        {
1240            chunksize=leftsize;
1241            rightstep=leftsize;
1242            leftstep=0;
1243            numsteps=m_right->getNumDPPSample();
1244            if (LS)
1245            {
1246               numsteps*=rightsize;
1247            }
1248        }
1249      }
1250    
1251      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1252      // Get the values of sub-expressions      // Get the values of sub-expressions
1253    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1254    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note      // calcBufss for why we can't put offset as the 2nd param above
1255      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1256      // the right child starts further along.      // the right child starts further along.
1257    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1258    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1259    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1260    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1261    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1262    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1263    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1264    
1265    
1266    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1267    switch(m_op)    switch(m_op)
1268    {    {
1269      case ADD:      case ADD:
1270      PROC_OP(plus<double>());          PROC_OP(NO_ARG,plus<double>());
1271      break;      break;
1272      case SUB:      case SUB:
1273      PROC_OP(minus<double>());      PROC_OP(NO_ARG,minus<double>());
1274      break;      break;
1275      case MUL:      case MUL:
1276      PROC_OP(multiplies<double>());      PROC_OP(NO_ARG,multiplies<double>());
1277      break;      break;
1278      case DIV:      case DIV:
1279      PROC_OP(divides<double>());      PROC_OP(NO_ARG,divides<double>());
1280      break;      break;
1281      case POW:      case POW:
1282      PROC_OP(::pow);         PROC_OP(double (double,double),::pow);
1283      break;      break;
1284      default:      default:
1285      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1286    }    }
1287    roffset=offset;      roffset=offset;
1288      return &v;
1289    }
1290    
1291    
1292    
1293    /*
1294      \brief Compute the value of the expression (tensor product) for the given sample.
1295      \return Vector which stores the value of the subexpression for the given sample.
1296      \param v A vector to store intermediate results.
1297      \param offset Index in v to begin storing results.
1298      \param sampleNo Sample number to evaluate.
1299      \param roffset (output parameter) the offset in the return vector where the result begins.
1300    
1301      The return value will be an existing vector so do not deallocate it.
1302      If the result is stored in v it should be stored at the offset given.
1303      Everything from offset to the end of v should be considered available for this method to use.
1304    */
1305    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1306    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1307    // unlike the other resolve helpers, we must treat these datapoints separately.
1308    DataTypes::ValueType*
1309    DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1310    {
1311    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1312    
1313      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1314        // first work out which of the children are expanded
1315      bool leftExp=(m_left->m_readytype=='E');
1316      bool rightExp=(m_right->m_readytype=='E');
1317      int steps=getNumDPPSample();
1318    /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1319      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1320      int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1321      int rightStep=(rightExp?m_right->getNoValues() : 0);
1322    
1323      int resultStep=getNoValues();
1324        // Get the values of sub-expressions (leave a gap of one sample for the result).
1325      int gap=offset+m_samplesize;  
1326    
1327    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1328    
1329      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1330      gap+=m_left->getMaxSampleSize();
1331    
1332    
1333    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1334    
1335    
1336      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1337    
1338    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1339    cout << getNoValues() << endl;)
1340    LAZYDEBUG(cerr << "Result of left=";)
1341    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1342    
1343    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1344    {
1345    cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1346    if (i%4==0) cout << endl;
1347    })
1348    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1349    LAZYDEBUG(
1350    for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1351    {
1352    cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1353    if (i%4==0) cout << endl;
1354    }
1355    cerr << endl;
1356    )
1357    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1358    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1359    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1360    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1361    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1362    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1363    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1364    
1365      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1366      switch(m_op)
1367      {
1368        case PROD:
1369        for (int i=0;i<steps;++i,resultp+=resultStep)
1370        {
1371    
1372    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1373    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1374    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1375    
1376              const double *ptr_0 = &((*left)[lroffset]);
1377              const double *ptr_1 = &((*right)[rroffset]);
1378    
1379    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1380    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1381    
1382              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1383    
1384    LAZYDEBUG(cout << "Results--\n";
1385    {
1386      DataVector dv(getNoValues());
1387    for (int z=0;z<getNoValues();++z)
1388    {
1389      cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1390      if (z%4==0) cout << endl;
1391      dv[z]=resultp[z];
1392    }
1393    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1394    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1395    }
1396    )
1397          lroffset+=leftStep;
1398          rroffset+=rightStep;
1399        }
1400        break;
1401        default:
1402        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1403      }
1404      roffset=offset;
1405    return &v;    return &v;
1406  }  }
1407    
# Line 715  cout << "Resolve binary: " << toString() Line 1426  cout << "Resolve binary: " << toString()
1426  const DataTypes::ValueType*  const DataTypes::ValueType*
1427  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1428  {  {
1429  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1430      // collapse so we have a 'E' node or an IDENTITY for some other type      // collapse so we have a 'E' node or an IDENTITY for some other type
1431    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1432    {    {
# Line 723  cout << "Resolve sample " << toString() Line 1434  cout << "Resolve sample " << toString()
1434    }    }
1435    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
1436    {    {
1437      const ValueType& vec=m_id->getVector();      const ValueType& vec=m_id->getVectorRO();
1438      if (m_readytype=='C')      if (m_readytype=='C')
1439      {      {
1440      roffset=0;      roffset=0;
1441    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1442      return &(vec);      return &(vec);
1443      }      }
1444      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1445    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1446      return &(vec);      return &(vec);
1447    }    }
1448    if (m_readytype!='E')    if (m_readytype!='E')
# Line 738  cout << "Resolve sample " << toString() Line 1451  cout << "Resolve sample " << toString()
1451    }    }
1452    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1453    {    {
1454    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY:
1455      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1456    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1457      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1458      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1459      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1460    default:    default:
1461      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1462    }    }
1463    
1464    }
1465    
1466    const DataTypes::ValueType*
1467    DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
1468    {
1469    #ifdef _OPENMP
1470        int tid=omp_get_thread_num();
1471    #else
1472        int tid=0;
1473    #endif
1474        return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
1475    }
1476    
1477    
1478    // This needs to do the work of the idenity constructor
1479    void
1480    DataLazy::resolveToIdentity()
1481    {
1482       if (m_op==IDENTITY)
1483        return;
1484       DataReady_ptr p=resolve();
1485       makeIdentity(p);
1486  }  }
1487    
1488    void DataLazy::makeIdentity(const DataReady_ptr& p)
1489    {
1490       m_op=IDENTITY;
1491       m_axis_offset=0;
1492       m_transpose=0;
1493       m_SL=m_SM=m_SR=0;
1494       m_children=m_height=0;
1495       m_id=p;
1496       if(p->isConstant()) {m_readytype='C';}
1497       else if(p->isExpanded()) {m_readytype='E';}
1498       else if (p->isTagged()) {m_readytype='T';}
1499       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1500       m_buffsRequired=1;
1501       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1502       m_maxsamplesize=m_samplesize;
1503       m_left.reset();
1504       m_right.reset();
1505    }
1506    
1507  // To simplify the memory management, all threads operate on one large vector, rather than one each.  // To simplify the memory management, all threads operate on one large vector, rather than one each.
1508  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
# Line 752  DataReady_ptr Line 1510  DataReady_ptr
1510  DataLazy::resolve()  DataLazy::resolve()
1511  {  {
1512    
1513  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1514  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
   
1515    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1516    {    {
1517      collapse();      collapse();
# Line 764  cout << "Buffers=" << m_buffsRequired << Line 1521  cout << "Buffers=" << m_buffsRequired <<
1521      return m_id;      return m_id;
1522    }    }
1523      // 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'
1524    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough    size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1525      // storage to evaluate its expression      // storage to evaluate its expression
1526    int numthreads=1;    int numthreads=1;
1527  #ifdef _OPENMP  #ifdef _OPENMP
1528    numthreads=getNumberOfThreads();    numthreads=omp_get_max_threads();
   int threadnum=0;  
1529  #endif  #endif
1530    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1531  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1532    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1533    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
1534    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1535    int sample;    int sample;
1536    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
1537    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1538    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1539    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1540    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1541      #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1542    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1543    {    {
1544  cout << "################################# " << sample << endl;        if (sample==0)  {ENABLEDEBUG}
1545    
1546    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1547    LAZYDEBUG(cout << "################################# " << sample << endl;)
1548  #ifdef _OPENMP  #ifdef _OPENMP
1549      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1550  #else  #else
1551      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
1552  #endif  #endif
1553  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1554    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1555      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1556  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1557      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
1558      {      {
1559    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1560      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1561      }      }
1562  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1563    LAZYDEBUG(cerr << "*********************************" << endl;)
1564        DISABLEDEBUG
1565    }    }
1566    return resptr;    return resptr;
1567  }  }
# Line 815  DataLazy::toString() const Line 1579  DataLazy::toString() const
1579  void  void
1580  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1581  {  {
1582    //    oss << "[" << m_children <<";"<<m_height <<"]";
1583    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1584    {    {
1585    case G_IDENTITY:    case G_IDENTITY:
# Line 844  DataLazy::intoString(ostringstream& oss) Line 1609  DataLazy::intoString(ostringstream& oss)
1609      oss << ')';      oss << ')';
1610      break;      break;
1611    case G_UNARY:    case G_UNARY:
1612      case G_UNARY_P:
1613      case G_NP1OUT:
1614      case G_NP1OUT_P:
1615      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1616      m_left->intoString(oss);      m_left->intoString(oss);
1617      oss << ')';      oss << ')';
1618      break;      break;
1619      case G_TENSORPROD:
1620        oss << opToString(m_op) << '(';
1621        m_left->intoString(oss);
1622        oss << ", ";
1623        m_right->intoString(oss);
1624        oss << ')';
1625        break;
1626    default:    default:
1627      oss << "UNKNOWN";      oss << "UNKNOWN";
1628    }    }
# Line 861  DataLazy::deepCopy() Line 1636  DataLazy::deepCopy()
1636    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1637    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1638    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1639      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1640      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1641    default:    default:
1642      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1643    }    }
1644  }  }
1645    
1646    
1647    // There is no single, natural interpretation of getLength on DataLazy.
1648    // Instances of DataReady can look at the size of their vectors.
1649    // For lazy though, it could be the size the data would be if it were resolved;
1650    // or it could be some function of the lengths of the DataReady instances which
1651    // form part of the expression.
1652    // Rather than have people making assumptions, I have disabled the method.
1653  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1654  DataLazy::getLength() const  DataLazy::getLength() const
1655  {  {
1656    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1657  }  }
1658    
1659    
# Line 936  DataLazy::getPointOffset(int sampleNo, Line 1719  DataLazy::getPointOffset(int sampleNo,
1719    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1720  }  }
1721    
1722  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
1723  // 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.  
1724  void  void
1725  DataLazy::setToZero()  DataLazy::setToZero()
1726  {  {
1727    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
1728    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1729    m_op=IDENTITY;  //   m_op=IDENTITY;
1730    m_right.reset();    //   m_right.reset();  
1731    m_left.reset();  //   m_left.reset();
1732    m_readytype='C';  //   m_readytype='C';
1733    m_buffsRequired=1;  //   m_buffsRequired=1;
1734    
1735      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1736    
1737    }
1738    
1739    bool
1740    DataLazy::actsExpanded() const
1741    {
1742        return (m_readytype=='E');
1743  }  }
1744    
1745  }   // end namespace  }   // end namespace

Legend:
Removed from v.1950  
changed lines
  Added in v.2271

  ViewVC Help
Powered by ViewVC 1.1.26