/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataLazy.cpp

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1898 by jfenwick, Mon Oct 20 01:20:18 2008 UTC trunk/escript/src/DataLazy.cpp revision 2037 by jfenwick, Thu Nov 13 06:17:12 2008 UTC
# 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    /*
32    How does DataLazy work?
33    ~~~~~~~~~~~~~~~~~~~~~~~
34    
35    Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36    denoted left and right.
37    
38    A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39    This means that all "internal" nodes in the structure are instances of DataLazy.
40    
41    Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42    Note that IDENITY is not considered a unary operation.
43    
44    I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45    It must however form a DAG (directed acyclic graph).
46    I will refer to individual DataLazy objects with the structure as nodes.
47    
48    Each node also stores:
49    - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50        evaluated.
51    - m_length ~ how many values would be stored in the answer if the expression was evaluated.
52    - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
53        evaluate the expression.
54    - m_samplesize ~ the number of doubles stored in a sample.
55    
56    When a new node is created, the above values are computed based on the values in the child nodes.
57    Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
58    
59    The resolve method, which produces a DataReady from a DataLazy, does the following:
60    1) Create a DataReady to hold the new result.
61    2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
62    3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
63    
64    (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
65    
66    resolveSample returns a Vector* and an offset within that vector where the result is stored.
67    Normally, this would be v, but for identity nodes their internal vector is returned instead.
68    
69    The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
70    
71    For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
72    The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
73    
74    To add a new operator you need to do the following (plus anything I might have forgotten):
75    1) Add to the ES_optype.
76    2) determine what opgroup your operation belongs to (X)
77    3) add a string for the op to the end of ES_opstrings
78    4) increase ES_opcount
79    5) add an entry (X) to opgroups
80    6) add an entry to the switch in collapseToReady
81    7) add an entry to resolveX
82    */
83    
84    
85  using namespace std;  using namespace std;
86  using namespace boost;  using namespace boost;
# Line 33  using namespace boost; Line 88  using namespace boost;
88  namespace escript  namespace escript
89  {  {
90    
 const std::string&  
 opToString(ES_optype op);  
   
91  namespace  namespace
92  {  {
93    
   
   
94  enum ES_opgroup  enum ES_opgroup
95  {  {
96     G_UNKNOWN,     G_UNKNOWN,
97     G_IDENTITY,     G_IDENTITY,
98     G_BINARY,     G_BINARY,        // pointwise operations with two arguments
99     G_UNARY     G_UNARY,     // pointwise operations with one argument
100       G_NP1OUT     // non-pointwise op with one output
101  };  };
102    
103    
104    
105    
106  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107                "sin","cos","tan",
108              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
109              "asinh","acosh","atanh",              "asinh","acosh","atanh",
110              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
111              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0",
112  int ES_opcount=32;              "symmetric","nonsymmetric"};
113  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9  int ES_opcount=35;
114              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
115              G_UNARY,G_UNARY,G_UNARY,                    // 19              G_UNARY,G_UNARY,G_UNARY, //10
116              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
117              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY,                    // 20
118                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
119                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,            // 33
120                G_NP1OUT,G_NP1OUT};
121  inline  inline
122  ES_opgroup  ES_opgroup
123  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 79  resultFS(DataAbstract_ptr left, DataAbst Line 134  resultFS(DataAbstract_ptr left, DataAbst
134      // 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
135      // programming error exception.      // programming error exception.
136    
137      FunctionSpace l=left->getFunctionSpace();
138      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
139      {    if (l!=r)
140          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
141      }      if (r.probeInterpolation(l))
142      return left->getFunctionSpace();      {
143        return l;
144        }
145        if (l.probeInterpolation(r))
146        {
147        return r;
148        }
149        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
150      }
151      return l;
152  }  }
153    
154  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
# Line 93  resultShape(DataAbstract_ptr left, DataA Line 157  resultShape(DataAbstract_ptr left, DataA
157  {  {
158      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
159      {      {
160          throw DataException("Shapes not the same - shapes must match for lazy data.");        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
161          {
162            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
163          }
164          if (left->getRank()==0)   // we need to allow scalar * anything
165          {
166            return right->getShape();
167          }
168          if (right->getRank()==0)
169          {
170            return left->getShape();
171          }
172          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
173      }      }
174      return left->getShape();      return left->getShape();
175  }  }
176    
177    // determine the number of points in the result of "left op right"
178  size_t  size_t
179  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
180  {  {
# Line 105  resultLength(DataAbstract_ptr left, Data Line 182  resultLength(DataAbstract_ptr left, Data
182     {     {
183     case G_BINARY: return left->getLength();     case G_BINARY: return left->getLength();
184     case G_UNARY: return left->getLength();     case G_UNARY: return left->getLength();
185       case G_NP1OUT: return left->getLength();
186     default:     default:
187      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
188     }     }
189  }  }
190    
191    // determine the number of samples requires to evaluate an expression combining left and right
192    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
193  int  int
194  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
195  {  {
# Line 118  calcBuffs(const DataLazy_ptr& left, cons Line 198  calcBuffs(const DataLazy_ptr& left, cons
198     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
199     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
200     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY: return max(left->getBuffsRequired(),1);
201       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
202     default:     default:
203      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
204     }     }
205  }  }
206    
207    
208  }   // end anonymous namespace  }   // end anonymous namespace
209    
210    
211    
212    // Return a string representing the operation
213  const std::string&  const std::string&
214  opToString(ES_optype op)  opToString(ES_optype op)
215  {  {
# Line 143  DataLazy::DataLazy(DataAbstract_ptr p) Line 227  DataLazy::DataLazy(DataAbstract_ptr p)
227  {  {
228     if (p->isLazy())     if (p->isLazy())
229     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
230      // I don't want identity of Lazy.      // I don't want identity of Lazy.
231      // Question: Why would that be so bad?      // Question: Why would that be so bad?
232      // 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 163  DataLazy::DataLazy(DataAbstract_ptr p) Line 246  DataLazy::DataLazy(DataAbstract_ptr p)
246  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
247  }  }
248    
249    
250    
251    
252  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
253      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
254      m_op(op)      m_op(op)
255  {  {
256     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
257     {     {
258      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.");
259     }     }
# Line 183  DataLazy::DataLazy(DataAbstract_ptr left Line 269  DataLazy::DataLazy(DataAbstract_ptr left
269     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
270     m_length=left->getLength();     m_length=left->getLength();
271     m_left=lleft;     m_left=lleft;
272     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
273     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
274  }  }
275    
276    
277  // DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // In this constructor we need to consider interpolation
 //  : parent(resultFS(left,right,op), resultShape(left,right,op)),  
 //  m_left(left),  
 //  m_right(right),  
 //  m_op(op)  
 // {  
 //    if (getOpgroup(op)!=G_BINARY)  
 //    {  
 //  throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");  
 //    }  
 //    m_length=resultLength(m_left,m_right,m_op);  
 //    m_samplesize=getNumDPPSample()*getNoValues();  
 //    m_buffsRequired=calcBuffs(m_left, m_right, m_op);  
 // cout << "(2)Lazy created with " << m_samplesize << endl;  
 // }  
   
278  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
279      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
280      m_op(op)      m_op(op)
281  {  {
282     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_BINARY))
283     {     {
284    cout << opToString(op) << endl;
285      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.");
286     }     }
287     if (left->isLazy())  
288       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
289       {
290        FunctionSpace fs=getFunctionSpace();
291        Data ltemp(left);
292        Data tmp(ltemp,fs);
293        left=tmp.borrowDataPtr();
294       }
295       if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated
296       {
297        Data tmp(Data(right),getFunctionSpace());
298        right=tmp.borrowDataPtr();
299       }
300       left->operandCheck(*right);
301    
302       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
303     {     {
304      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
305     }     }
# Line 243  DataLazy::DataLazy(DataAbstract_ptr left Line 330  DataLazy::DataLazy(DataAbstract_ptr left
330      m_readytype='C';      m_readytype='C';
331     }     }
332     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
333     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
334     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
335  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
336  }  }
# Line 261  DataLazy::getBuffsRequired() const Line 348  DataLazy::getBuffsRequired() const
348  }  }
349    
350    
351    /*
352      \brief Evaluates the expression using methods on Data.
353      This does the work for the collapse method.
354      For reasons of efficiency do not call this method on DataExpanded nodes.
355    */
356  DataReady_ptr  DataReady_ptr
357  DataLazy::collapseToReady()  DataLazy::collapseToReady()
358  {  {
# Line 374  DataLazy::collapseToReady() Line 466  DataLazy::collapseToReady()
466      case LEZ:      case LEZ:
467      result=left.whereNonPositive();      result=left.whereNonPositive();
468      break;      break;
469        case SYM:
470        result=left.symmetric();
471        break;
472        case NSYM:
473        result=left.nonsymmetric();
474        break;
475      default:      default:
476      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)+".");
477    }    }
478    return result.borrowReadyPtr();    return result.borrowReadyPtr();
479  }  }
480    
481    /*
482       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
483       This method uses the original methods on the Data class to evaluate the expressions.
484       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
485       the purpose of using DataLazy in the first place).
486    */
487  void  void
488  DataLazy::collapse()  DataLazy::collapse()
489  {  {
# Line 395  DataLazy::collapse() Line 499  DataLazy::collapse()
499    m_op=IDENTITY;    m_op=IDENTITY;
500  }  }
501    
502    /*
503      \brief Compute the value of the expression (unary operation) for the given sample.
504      \return Vector which stores the value of the subexpression for the given sample.
505      \param v A vector to store intermediate results.
506      \param offset Index in v to begin storing results.
507      \param sampleNo Sample number to evaluate.
508      \param roffset (output parameter) the offset in the return vector where the result begins.
509    
510      The return value will be an existing vector so do not deallocate it.
511      If the result is stored in v it should be stored at the offset given.
512      Everything from offset to the end of v should be considered available for this method to use.
513    */
514  DataTypes::ValueType*  DataTypes::ValueType*
515  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
516  {  {
# Line 412  DataLazy::resolveUnary(ValueType& v, siz Line 528  DataLazy::resolveUnary(ValueType& v, siz
528    switch (m_op)    switch (m_op)
529    {    {
530      case SIN:        case SIN:  
531      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
532      break;      break;
533      case COS:      case COS:
534      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
535      break;      break;
536      case TAN:      case TAN:
537      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
538      break;      break;
539      case ASIN:      case ASIN:
540      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
541      break;      break;
542      case ACOS:      case ACOS:
543      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
544      break;      break;
545      case ATAN:      case ATAN:
546      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
547      break;      break;
548      case SINH:      case SINH:
549      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
550      break;      break;
551      case COSH:      case COSH:
552      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
553      break;      break;
554      case TANH:      case TANH:
555      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
556      break;      break;
557      case ERF:      case ERF:
558  #ifdef _WIN32  #ifdef _WIN32
# Line 467  DataLazy::resolveUnary(ValueType& v, siz Line 583  DataLazy::resolveUnary(ValueType& v, siz
583  #endif    #endif  
584      break;      break;
585      case LOG10:      case LOG10:
586      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
587      break;      break;
588      case LOG:      case LOG:
589      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
590      break;      break;
591      case SIGN:      case SIGN:
592      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
593      break;      break;
594      case ABS:      case ABS:
595      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
596      break;      break;
597      case NEG:      case NEG:
598      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 487  DataLazy::resolveUnary(ValueType& v, siz Line 603  DataLazy::resolveUnary(ValueType& v, siz
603      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
604      break;      break;
605      case EXP:      case EXP:
606      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
607      break;      break;
608      case SQRT:      case SQRT:
609      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
610      break;      break;
611      case RECIP:      case RECIP:
612      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 515  DataLazy::resolveUnary(ValueType& v, siz Line 631  DataLazy::resolveUnary(ValueType& v, siz
631  }  }
632    
633    
634    /*
635      \brief Compute the value of the expression (unary operation) for the given sample.
636      \return Vector which stores the value of the subexpression for the given sample.
637      \param v A vector to store intermediate results.
638      \param offset Index in v to begin storing results.
639      \param sampleNo Sample number to evaluate.
640      \param roffset (output parameter) the offset in the return vector where the result begins.
641    
642      The return value will be an existing vector so do not deallocate it.
643      If the result is stored in v it should be stored at the offset given.
644      Everything from offset to the end of v should be considered available for this method to use.
645    */
646    DataTypes::ValueType*
647    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
648    {
649        // we assume that any collapsing has been done before we get here
650        // since we only have one argument we don't need to think about only
651        // processing single points.
652      if (m_readytype!='E')
653      {
654        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
655      }
656        // since we can't write the result over the input, we need a result offset further along
657      size_t subroffset=roffset+m_samplesize;
658      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
659      roffset=offset;
660      switch (m_op)
661      {
662        case SYM:
663        DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
664        break;
665        case NSYM:
666        DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
667        break;
668        default:
669        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
670      }
671      return &v;
672    }
673    
674    
675    
 // const double*  
 // DataLazy::resolveUnary(ValueType& v,int sampleNo,  size_t offset) const  
 // {  
 //  // we assume that any collapsing has been done before we get here  
 //  // since we only have one argument we don't need to think about only  
 //  // processing single points.  
 //   if (m_readytype!='E')  
 //   {  
 //     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");  
 //   }  
 //   const double* left=m_left->resolveSample(v,sampleNo,offset);  
 //   double* result=&(v[offset]);  
 //   switch (m_op)  
 //   {  
 //     case SIN:      
 //  tensor_unary_operation(m_samplesize, left, result, ::sin);  
 //  break;  
 //     case COS:  
 //  tensor_unary_operation(m_samplesize, left, result, ::cos);  
 //  break;  
 //     case TAN:  
 //  tensor_unary_operation(m_samplesize, left, result, ::tan);  
 //  break;  
 //     case ASIN:  
 //  tensor_unary_operation(m_samplesize, left, result, ::asin);  
 //  break;  
 //     case ACOS:  
 //  tensor_unary_operation(m_samplesize, left, result, ::acos);  
 //  break;  
 //     case ATAN:  
 //  tensor_unary_operation(m_samplesize, left, result, ::atan);  
 //  break;  
 //     case SINH:  
 //  tensor_unary_operation(m_samplesize, left, result, ::sinh);  
 //  break;  
 //     case COSH:  
 //  tensor_unary_operation(m_samplesize, left, result, ::cosh);  
 //  break;  
 //     case TANH:  
 //  tensor_unary_operation(m_samplesize, left, result, ::tanh);  
 //  break;  
 //     case ERF:  
 // #ifdef _WIN32  
 //  throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::erf);  
 //  break;  
 // #endif  
 //    case ASINH:  
 // #ifdef _WIN32  
 //  tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 // #endif    
 //  break;  
 //    case ACOSH:  
 // #ifdef _WIN32  
 //  tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 // #endif    
 //  break;  
 //    case ATANH:  
 // #ifdef _WIN32  
 //  tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 // #endif    
 //  break;  
 //     case LOG10:  
 //  tensor_unary_operation(m_samplesize, left, result, ::log10);  
 //  break;  
 //     case LOG:  
 //  tensor_unary_operation(m_samplesize, left, result, ::log);  
 //  break;  
 //     case SIGN:  
 //  tensor_unary_operation(m_samplesize, left, result, escript::fsign);  
 //  break;  
 //     case ABS:  
 //  tensor_unary_operation(m_samplesize, left, result, ::fabs);  
 //  break;  
 //     case NEG:  
 //  tensor_unary_operation(m_samplesize, left, result, negate<double>());  
 //  break;  
 //     case POS:  
 //  // it doesn't mean anything for delayed.  
 //  // it will just trigger a deep copy of the lazy object  
 //  throw DataException("Programmer error - POS not supported for lazy data.");  
 //  break;  
 //     case EXP:  
 //  tensor_unary_operation(m_samplesize, left, result, ::exp);  
 //  break;  
 //     case SQRT:  
 //  tensor_unary_operation(m_samplesize, left, result, ::sqrt);  
 //  break;  
 //     case RECIP:  
 //  tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));  
 //  break;  
 //     case GZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));  
 //  break;  
 //     case LZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));  
 //  break;  
 //     case GEZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));  
 //  break;  
 //     case LEZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));  
 //  break;  
 //  
 //     default:  
 //  throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");  
 //   }  
 //   return result;  
 // }  
676    
677  #define PROC_OP(X) \  #define PROC_OP(TYPE,X)                               \
678      for (int i=0;i<steps;++i,resultp+=getNoValues()) \      for (int i=0;i<steps;++i,resultp+=resultStep) \
679      { \      { \
680  cout << "Step#" << i << " chunk=" << chunksize << endl; \         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
681  cout << left[0] << left[1] << left[2] << endl; \         lroffset+=leftStep; \
682  cout << right[0] << right[1] << right[2] << endl; \         rroffset+=rightStep; \
        tensor_binary_operation(chunksize, left, right, resultp, X); \  
        left+=leftStep; \  
        right+=rightStep; \  
 cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \  
683      }      }
684    
685    /*
686      \brief Compute the value of the expression (binary operation) for the given sample.
687      \return Vector which stores the value of the subexpression for the given sample.
688      \param v A vector to store intermediate results.
689      \param offset Index in v to begin storing results.
690      \param sampleNo Sample number to evaluate.
691      \param roffset (output parameter) the offset in the return vector where the result begins.
692    
693      The return value will be an existing vector so do not deallocate it.
694      If the result is stored in v it should be stored at the offset given.
695      Everything from offset to the end of v should be considered available for this method to use.
696    */
697    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
698    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
699    // If both children are expanded, then we can process them in a single operation (we treat
700    // the whole sample as one big datapoint.
701    // If one of the children is not expanded, then we need to treat each point in the sample
702    // individually.
703    // There is an additional complication when scalar operations are considered.
704    // For example, 2+Vector.
705    // In this case each double within the point is treated individually
706  DataTypes::ValueType*  DataTypes::ValueType*
707  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
708  {  {
     // again we assume that all collapsing has already been done  
     // so we have at least one expanded child.  
     // however, we could still have one of the children being not expanded.  
   
709  cout << "Resolve binary: " << toString() << endl;  cout << "Resolve binary: " << toString() << endl;
710    
711    size_t lroffset=0, rroffset=0;    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
712        // first work out which of the children are expanded
713    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
714    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
715    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
716    int steps=(bigloops?1:getNumDPPSample());    int steps=(bigloops?1:getNumDPPSample());
717    size_t chunksize=(bigloops? m_samplesize : getNoValues());    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
718    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
719    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);    {
720        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
721        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
722        chunksize=1;    // for scalar
723      }    
724      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
725      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
726      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
727        // Get the values of sub-expressions
728    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
729    const ValueType* right=m_right->resolveSample(v,offset,sampleNo,rroffset);        const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
730      // now we need to know which args are expanded      // the right child starts further along.
731  cout << "left=" << left << " right=" << right << endl;    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
 cout << "(Length) l=" << left->size() << " r=" << right->size() << " res=" << v.size() << endl;  
   double* resultp=&(v[offset]);  
732    switch(m_op)    switch(m_op)
733    {    {
734      case ADD:      case ADD:
735      for (int i=0;i<steps;++i,resultp+=getNoValues())          PROC_OP(NO_ARG,plus<double>());
736      {      break;
737  cerr << "Step#" << i << " chunk=" << chunksize << endl;      case SUB:
738  cerr << left << "[" << lroffset << "] " << right << "[" << rroffset << "]" << endl;      PROC_OP(NO_ARG,minus<double>());
739         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, plus<double>());      break;
740         lroffset+=leftStep;      case MUL:
741         rroffset+=rightStep;      PROC_OP(NO_ARG,multiplies<double>());
742  cerr << "left=" << lroffset << " right=" << rroffset << endl;      break;
743      }      case DIV:
744        PROC_OP(NO_ARG,divides<double>());
745        break;
746        case POW:
747           PROC_OP(double (double,double),::pow);
748      break;      break;
 // need to fill in the rest  
749      default:      default:
750      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
751    }    }
752    roffset=offset;    roffset=offset;  
753    return &v;    return &v;
754  }  }
755    
756    
757    
758  // #define PROC_OP(X) \  /*
759  //  for (int i=0;i<steps;++i,resultp+=getNoValues()) \    \brief Compute the value of the expression for the given sample.
760  //  { \    \return Vector which stores the value of the subexpression for the given sample.
761  // cout << "Step#" << i << " chunk=" << chunksize << endl; \    \param v A vector to store intermediate results.
762  // cout << left[0] << left[1] << left[2] << endl; \    \param offset Index in v to begin storing results.
763  // cout << right[0] << right[1] << right[2] << endl; \    \param sampleNo Sample number to evaluate.
764  //     tensor_binary_operation(chunksize, left, right, resultp, X); \    \param roffset (output parameter) the offset in the return vector where the result begins.
 //     left+=leftStep; \  
 //     right+=rightStep; \  
 // cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \  
 //  }  
 //  
 // const double*  
 // DataLazy::resolveBinary(ValueType& v,int sampleNo,  size_t offset) const  
 // {  
 //  // again we assume that all collapsing has already been done  
 //  // so we have at least one expanded child.  
 //  // however, we could still have one of the children being not expanded.  
 //  
 // cout << "Resolve binary: " << toString() << endl;  
 //  
 //   const double* left=m_left->resolveSample(v,sampleNo,offset);  
 // // cout << "Done Left " << /*left[0] << left[1] << left[2] << */endl;  
 //   const double* right=m_right->resolveSample(v,sampleNo,offset);  
 // // cout << "Done Right"  << /*right[0] << right[1] << right[2] <<*/ endl;  
 //      // now we need to know which args are expanded  
 //   bool leftExp=(m_left->m_readytype=='E');  
 //   bool rightExp=(m_right->m_readytype=='E');  
 //   bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp));  // is processing in single step  
 //   int steps=(bigloops?1:getNumSamples());  
 //   size_t chunksize=(bigloops? m_samplesize : getNoValues());  
 //   int leftStep=((leftExp && !rightExp)? getNoValues() : 0);  
 //   int rightStep=((rightExp && !leftExp)? getNoValues() : 0);  
 // cout << "left=" << left << " right=" << right << endl;  
 //   double* result=&(v[offset]);  
 //   double* resultp=result;  
 //   switch(m_op)  
 //   {  
 //     case ADD:  
 //  for (int i=0;i<steps;++i,resultp+=getNoValues())  
 //  {  
 // cout << "Step#" << i << " chunk=" << chunksize << endl; \  
 // // cout << left[0] << left[1] << left[2] << endl;  
 // // cout << right[0] << right[1] << right[2] << endl;  
 //     tensor_binary_operation(chunksize, left, right, resultp, plus<double>());  
 // cout << "left=" << left << " right=" << right << " resp=" << resultp << endl;  
 //     left+=leftStep;  
 //     right+=rightStep;  
 // cout << "left=" << left << " right=" << right << endl;  
 // // cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl;  
 //  }  
 //  break;  
 // // need to fill in the rest  
 //     default:  
 //  throw DataException("Programmer error - resolveBinay can not resolve operator "+opToString(m_op)+".");  
 //   }  
 // // cout << "About to return "  << result[0] << result[1] << result[2] << endl;;  
 //   return result;  
 // }  
   
 // // the vector and the offset are a place where the method could write its data if it wishes  
 // // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
 // // Hence the return value to indicate where the data is actually stored.  
 // // Regardless, the storage should be assumed to be used, even if it isn't.  
 // const double*  
 // DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset )  
 // {  
 // cout << "Resolve sample " << toString() << endl;  
 //  // collapse so we have a 'E' node or an IDENTITY for some other type  
 //   if (m_readytype!='E' && m_op!=IDENTITY)  
 //   {  
 //  collapse();  
 //   }  
 //   if (m_op==IDENTITY)      
 //   {  
 //     const ValueType& vec=m_id->getVector();  
 //     if (m_readytype=='C')  
 //     {  
 //  return &(vec[0]);  
 //     }  
 //     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
 //   }  
 //   if (m_readytype!='E')  
 //   {  
 //     throw DataException("Programmer Error - Collapse did not produce an expanded node.");  
 //   }  
 //   switch (getOpgroup(m_op))  
 //   {  
 //   case G_UNARY: return resolveUnary(v,sampleNo,offset);  
 //   case G_BINARY: return resolveBinary(v,sampleNo,offset);  
 //   default:  
 //     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
 //   }  
 // }  
   
   
765    
766      The return value will be an existing vector so do not deallocate it.
767    */
768  // the vector and the offset are a place where the method could write its data if it wishes  // the vector and the offset are a place where the method could write its data if it wishes
769  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  // it is not obligated to do so. For example, if it has its own storage already, it can use that.
770  // Hence the return value to indicate where the data is actually stored.  // Hence the return value to indicate where the data is actually stored.
# Line 830  cout << "Resolve sample " << toString() Line 805  cout << "Resolve sample " << toString()
805  }  }
806    
807    
808    // To simplify the memory management, all threads operate on one large vector, rather than one each.
809    // Each sample is evaluated independently and copied into the result DataExpanded.
 // This version uses double* trying again with vectors  
 // DataReady_ptr  
 // DataLazy::resolve()  
 // {  
 //  
 // cout << "Sample size=" << m_samplesize << endl;  
 // cout << "Buffers=" << m_buffsRequired << endl;  
 //  
 //   if (m_readytype!='E')  
 //   {  
 //     collapse();  
 //   }  
 //   if (m_op==IDENTITY)  
 //   {  
 //     return m_id;  
 //   }  
 //      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'  
 //   size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);  
 //   int numthreads=1;  
 // #ifdef _OPENMP  
 //   numthreads=getNumberOfThreads();  
 //   int threadnum=0;  
 // #endif  
 //   ValueType v(numthreads*threadbuffersize);    
 // cout << "Buffer created with size=" << v.size() << endl;  
 //   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));  
 //   ValueType& resvec=result->getVector();  
 //   DataReady_ptr resptr=DataReady_ptr(result);  
 //   int sample;  
 //   int resoffset;  
 //   int totalsamples=getNumSamples();  
 //   const double* res=0;  
 //   #pragma omp parallel for private(sample,resoffset,threadnum,res) schedule(static)  
 //   for (sample=0;sample<totalsamples;++sample)  
 //   {  
 // cout << "################################# " << sample << endl;  
 // #ifdef _OPENMP  
 //     res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());  
 // #else  
 //     res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.  
 // #endif  
 // cerr << "-------------------------------- " << endl;  
 //     resoffset=result->getPointOffset(sample,0);  
 // cerr << "offset=" << resoffset << endl;  
 //     for (unsigned int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector  
 //     {  
 //  resvec[resoffset]=res[i];  
 //     }  
 // cerr << "*********************************" << endl;  
 //   }  
 //   return resptr;  
 // }  
   
   
810  DataReady_ptr  DataReady_ptr
811  DataLazy::resolve()  DataLazy::resolve()
812  {  {
# Line 893  DataLazy::resolve() Line 814  DataLazy::resolve()
814  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
815  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
816    
817    if (m_readytype!='E')    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
818    {    {
819      collapse();      collapse();
820    }    }
821    if (m_op==IDENTITY)    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
822    {    {
823      return m_id;      return m_id;
824    }    }
825      // 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'
826    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough
827        // storage to evaluate its expression
828    int numthreads=1;    int numthreads=1;
829  #ifdef _OPENMP  #ifdef _OPENMP
830    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
# Line 916  cout << "Buffer created with size=" << v Line 838  cout << "Buffer created with size=" << v
838    int sample;    int sample;
839    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
840    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
841    const ValueType* res=0;    const ValueType* res=0;   // Vector storing the answer
842    size_t resoffset=0;    size_t resoffset=0;       // where in the vector to find the answer
843    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
844    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
845    {    {
# Line 925  cout << "############################### Line 847  cout << "###############################
847  #ifdef _OPENMP  #ifdef _OPENMP
848      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
849  #else  #else
850      res=resolveSample(v,0,sample,resoffset);   // this 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.
851  #endif  #endif
852  cerr << "-------------------------------- " << endl;  cerr << "-------------------------------- " << endl;
853      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
# Line 948  DataLazy::toString() const Line 870  DataLazy::toString() const
870    return oss.str();    return oss.str();
871  }  }
872    
873    
874  void  void
875  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
876  {  {
# Line 980  DataLazy::intoString(ostringstream& oss) Line 903  DataLazy::intoString(ostringstream& oss)
903      oss << ')';      oss << ')';
904      break;      break;
905    case G_UNARY:    case G_UNARY:
906      case G_NP1OUT:
907      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
908      m_left->intoString(oss);      m_left->intoString(oss);
909      oss << ')';      oss << ')';
# Line 989  DataLazy::intoString(ostringstream& oss) Line 913  DataLazy::intoString(ostringstream& oss)
913    }    }
914  }  }
915    
 // Note that in this case, deepCopy does not make copies of the leaves.  
 // Hopefully copy on write (or whatever we end up using) will take care of this.  
916  DataAbstract*  DataAbstract*
917  DataLazy::deepCopy()  DataLazy::deepCopy()
918  {  {
919    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
920    {    {
921      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
922      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
923      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
924      default:
925        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
926    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
927  }  }
928    
929    
# Line 1015  DataLazy::getSlice(const DataTypes::Regi Line 940  DataLazy::getSlice(const DataTypes::Regi
940    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
941  }  }
942    
943    
944    // To do this we need to rely on our child nodes
945    DataTypes::ValueType::size_type
946    DataLazy::getPointOffset(int sampleNo,
947                     int dataPointNo)
948    {
949      if (m_op==IDENTITY)
950      {
951        return m_id->getPointOffset(sampleNo,dataPointNo);
952      }
953      if (m_readytype!='E')
954      {
955        collapse();
956        return m_id->getPointOffset(sampleNo,dataPointNo);
957      }
958      // at this point we do not have an identity node and the expression will be Expanded
959      // so we only need to know which child to ask
960      if (m_left->m_readytype=='E')
961      {
962        return m_left->getPointOffset(sampleNo,dataPointNo);
963      }
964      else
965      {
966        return m_right->getPointOffset(sampleNo,dataPointNo);
967      }
968    }
969    
970    // To do this we need to rely on our child nodes
971  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
972  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
973                   int dataPointNo) const                   int dataPointNo) const
974  {  {
975    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
976      {
977        return m_id->getPointOffset(sampleNo,dataPointNo);
978      }
979      if (m_readytype=='E')
980      {
981        // at this point we do not have an identity node and the expression will be Expanded
982        // so we only need to know which child to ask
983        if (m_left->m_readytype=='E')
984        {
985        return m_left->getPointOffset(sampleNo,dataPointNo);
986        }
987        else
988        {
989        return m_right->getPointOffset(sampleNo,dataPointNo);
990        }
991      }
992      if (m_readytype=='C')
993      {
994        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
995      }
996      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
997    }
998    
999    // It would seem that DataTagged will need to be treated differently since even after setting all tags
1000    // to zero, all the tags from all the DataTags would be in the result.
1001    // However since they all have the same value (0) whether they are there or not should not matter.
1002    // So I have decided that for all types this method will create a constant 0.
1003    // It can be promoted up as required.
1004    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1005    // but we can deal with that if it arrises.
1006    void
1007    DataLazy::setToZero()
1008    {
1009      DataTypes::ValueType v(getNoValues(),0);
1010      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1011      m_op=IDENTITY;
1012      m_right.reset();  
1013      m_left.reset();
1014      m_readytype='C';
1015      m_buffsRequired=1;
1016  }  }
1017    
1018  }   // end namespace  }   // end namespace

Legend:
Removed from v.1898  
changed lines
  Added in v.2037

  ViewVC Help
Powered by ViewVC 1.1.26