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

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

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1910 by jfenwick, Thu Oct 23 03:05:28 2008 UTC trunk/escript/src/DataLazy.cpp revision 2177 by jfenwick, Wed Dec 17 23:51:23 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    // #define LAZYDEBUG(X) if (privdebug){X;}
32    #define LAZYDEBUG(X)
33    namespace
34    {
35    bool privdebug=false;
36    
37    #define ENABLEDEBUG privdebug=true;
38    #define DISABLEDEBUG privdebug=false;
39    }
40    
41    #define SIZELIMIT
42    // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
43    
44    
45  /*  /*
46  How does DataLazy work?  How does DataLazy work?
# Line 47  I will refer to individual DataLazy obje Line 62  I will refer to individual DataLazy obje
62  Each node also stores:  Each node also stores:
63  - 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
64      evaluated.      evaluated.
 - m_length ~ how many values would be stored in the answer if the expression was evaluated.  
65  - 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
66      evaluate the expression.      evaluate the expression.
67  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
# Line 69  The convention that I use, is that the r Line 83  The convention that I use, is that the r
83    
84  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.
85  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.
86    
87    To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
88    1) Add to the ES_optype.
89    2) determine what opgroup your operation belongs to (X)
90    3) add a string for the op to the end of ES_opstrings
91    4) increase ES_opcount
92    5) add an entry (X) to opgroups
93    6) add an entry to the switch in collapseToReady
94    7) add an entry to resolveX
95  */  */
96    
97    
# Line 78  using namespace boost; Line 101  using namespace boost;
101  namespace escript  namespace escript
102  {  {
103    
 const std::string&  
 opToString(ES_optype op);  
   
104  namespace  namespace
105  {  {
106    
# Line 89  enum ES_opgroup Line 109  enum ES_opgroup
109     G_UNKNOWN,     G_UNKNOWN,
110     G_IDENTITY,     G_IDENTITY,
111     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
112     G_UNARY      // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
113       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
114       G_NP1OUT,        // non-pointwise op with one output
115       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
116       G_TENSORPROD     // general tensor product
117  };  };
118    
119    
# Line 100  string ES_opstrings[]={"UNKNOWN","IDENTI Line 124  string ES_opstrings[]={"UNKNOWN","IDENTI
124              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
125              "asinh","acosh","atanh",              "asinh","acosh","atanh",
126              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
127              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
128  int ES_opcount=32;              "symmetric","nonsymmetric",
129                "prod",
130                "transpose", "trace"};
131    int ES_opcount=40;
132  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,
133              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
134              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
135              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
136              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
137              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
138                G_NP1OUT,G_NP1OUT,
139                G_TENSORPROD,
140                G_NP1OUT_P, G_NP1OUT_P};
141  inline  inline
142  ES_opgroup  ES_opgroup
143  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 124  resultFS(DataAbstract_ptr left, DataAbst Line 154  resultFS(DataAbstract_ptr left, DataAbst
154      // 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
155      // programming error exception.      // programming error exception.
156    
157      FunctionSpace l=left->getFunctionSpace();
158      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
159      {    if (l!=r)
160          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
161      }      if (r.probeInterpolation(l))
162      return left->getFunctionSpace();      {
163        return l;
164        }
165        if (l.probeInterpolation(r))
166        {
167        return r;
168        }
169        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
170      }
171      return l;
172  }  }
173    
174  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
175    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
176  DataTypes::ShapeType  DataTypes::ShapeType
177  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
178  {  {
179      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
180      {      {
181        if (getOpgroup(op)!=G_BINARY)        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
182        {        {
183          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.");
184        }        }
# Line 155  resultShape(DataAbstract_ptr left, DataA Line 195  resultShape(DataAbstract_ptr left, DataA
195      return left->getShape();      return left->getShape();
196  }  }
197    
198  // determine the number of points in the result of "left op right"  // return the shape for "op left"
199  size_t  
200  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataTypes::ShapeType
201    resultShape(DataAbstract_ptr left, ES_optype op)
202  {  {
203     switch (getOpgroup(op))      switch(op)
204     {      {
205     case G_BINARY: return left->getLength();          case TRANS:
206     case G_UNARY: return left->getLength();          return left->getShape();
207     default:      break;
208      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      case TRACE:
209     }          return DataTypes::scalarShape;
210        break;
211            default:
212        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
213        }
214    }
215    
216    // determine the output shape for the general tensor product operation
217    // the additional parameters return information required later for the product
218    // the majority of this code is copy pasted from C_General_Tensor_Product
219    DataTypes::ShapeType
220    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
221    {
222        
223      // Get rank and shape of inputs
224      int rank0 = left->getRank();
225      int rank1 = right->getRank();
226      const DataTypes::ShapeType& shape0 = left->getShape();
227      const DataTypes::ShapeType& shape1 = right->getShape();
228    
229      // Prepare for the loops of the product and verify compatibility of shapes
230      int start0=0, start1=0;
231      if (transpose == 0)       {}
232      else if (transpose == 1)  { start0 = axis_offset; }
233      else if (transpose == 2)  { start1 = rank1-axis_offset; }
234      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
235    
236      if (rank0<axis_offset)
237      {
238        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
239      }
240    
241      // Adjust the shapes for transpose
242      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
243      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
244      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
245      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
246    
247      // Prepare for the loops of the product
248      SL=1, SM=1, SR=1;
249      for (int i=0; i<rank0-axis_offset; i++)   {
250        SL *= tmpShape0[i];
251      }
252      for (int i=rank0-axis_offset; i<rank0; i++)   {
253        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
254          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
255        }
256        SM *= tmpShape0[i];
257      }
258      for (int i=axis_offset; i<rank1; i++)     {
259        SR *= tmpShape1[i];
260      }
261    
262      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
263      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
264      {         // block to limit the scope of out_index
265         int out_index=0;
266         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
267         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
268      }
269    
270      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
271      {
272         ostringstream os;
273         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
274         throw DataException(os.str());
275      }
276    
277      return shape2;
278  }  }
279    
280  // 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
281    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
282    // The same goes for G_TENSORPROD
283    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefths.
284    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
285    // multiple times
286  int  int
287  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
288  {  {
289     switch(getOpgroup(op))     switch(getOpgroup(op))
290     {     {
291     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
292     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
293     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY:
294       case G_UNARY_P: return max(left->getBuffsRequired(),1);
295       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
296       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
297       case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
298     default:     default:
299      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
300     }     }
# Line 200  opToString(ES_optype op) Line 318  opToString(ES_optype op)
318    
319    
320  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
321      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
     m_op(IDENTITY)  
322  {  {
323     if (p->isLazy())     if (p->isLazy())
324     {     {
# Line 212  DataLazy::DataLazy(DataAbstract_ptr p) Line 329  DataLazy::DataLazy(DataAbstract_ptr p)
329     }     }
330     else     else
331     {     {
332      m_id=dynamic_pointer_cast<DataReady>(p);      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
333      if(p->isConstant()) {m_readytype='C';}      makeIdentity(dr);
     else if(p->isExpanded()) {m_readytype='E';}  
     else if (p->isTagged()) {m_readytype='T';}  
     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}  
334     }     }
335     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;  
336  }  }
337    
338    
# Line 229  cout << "(1)Lazy created with " << m_sam Line 340  cout << "(1)Lazy created with " << m_sam
340    
341  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
342      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
343      m_op(op)      m_op(op),
344        m_axis_offset(0),
345        m_transpose(0),
346        m_SL(0), m_SM(0), m_SR(0)
347  {  {
348     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
349     {     {
350      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.");
351     }     }
352    
353     DataLazy_ptr lleft;     DataLazy_ptr lleft;
354     if (!left->isLazy())     if (!left->isLazy())
355     {     {
# Line 245  DataLazy::DataLazy(DataAbstract_ptr left Line 360  DataLazy::DataLazy(DataAbstract_ptr left
360      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
361     }     }
362     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
    m_length=left->getLength();  
363     m_left=lleft;     m_left=lleft;
364     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
365     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
366       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
367       m_children=m_left->m_children+1;
368       m_height=m_left->m_height+1;
369       SIZELIMIT
370  }  }
371    
372    
373    // In this constructor we need to consider interpolation
374  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
375      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
376      m_op(op)      m_op(op),
377        m_SL(0), m_SM(0), m_SR(0)
378  {  {
379     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_BINARY))
380     {     {
381      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.");
382     }     }
383    
384       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
385       {
386        FunctionSpace fs=getFunctionSpace();
387        Data ltemp(left);
388        Data tmp(ltemp,fs);
389        left=tmp.borrowDataPtr();
390       }
391       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
392       {
393        Data tmp(Data(right),getFunctionSpace());
394        right=tmp.borrowDataPtr();
395       }
396       left->operandCheck(*right);
397    
398     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
399     {     {
400      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
# Line 290  DataLazy::DataLazy(DataAbstract_ptr left Line 425  DataLazy::DataLazy(DataAbstract_ptr left
425     {     {
426      m_readytype='C';      m_readytype='C';
427     }     }
428     m_length=resultLength(m_left,m_right,m_op);     m_samplesize=getNumDPPSample()*getNoValues();
429     m_samplesize=getNumDPPSample()*getNoValues();         m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
430     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
431  cout << "(3)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
432       m_height=max(m_left->m_height,m_right->m_height)+1;
433       SIZELIMIT
434    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
435  }  }
436    
437    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
438        : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
439        m_op(op),
440        m_axis_offset(axis_offset),
441        m_transpose(transpose)
442    {
443       if ((getOpgroup(op)!=G_TENSORPROD))
444       {
445        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
446       }
447       if ((transpose>2) || (transpose<0))
448       {
449        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
450       }
451       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
452       {
453        FunctionSpace fs=getFunctionSpace();
454        Data ltemp(left);
455        Data tmp(ltemp,fs);
456        left=tmp.borrowDataPtr();
457       }
458       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
459       {
460        Data tmp(Data(right),getFunctionSpace());
461        right=tmp.borrowDataPtr();
462       }
463       left->operandCheck(*right);
464    
465       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
466       {
467        m_left=dynamic_pointer_cast<DataLazy>(left);
468       }
469       else
470       {
471        m_left=DataLazy_ptr(new DataLazy(left));
472       }
473       if (right->isLazy())
474       {
475        m_right=dynamic_pointer_cast<DataLazy>(right);
476       }
477       else
478       {
479        m_right=DataLazy_ptr(new DataLazy(right));
480       }
481       char lt=m_left->m_readytype;
482       char rt=m_right->m_readytype;
483       if (lt=='E' || rt=='E')
484       {
485        m_readytype='E';
486       }
487       else if (lt=='T' || rt=='T')
488       {
489        m_readytype='T';
490       }
491       else
492       {
493        m_readytype='C';
494       }
495       m_samplesize=getNumDPPSample()*getNoValues();
496       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
497       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
498       m_children=m_left->m_children+m_right->m_children+2;
499       m_height=max(m_left->m_height,m_right->m_height)+1;
500       SIZELIMIT
501    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
502    }
503    
504    
505    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
506        : parent(left->getFunctionSpace(), resultShape(left,op)),
507        m_op(op),
508        m_axis_offset(axis_offset),
509        m_transpose(0),
510        m_tol(0)
511    {
512       if ((getOpgroup(op)!=G_NP1OUT_P))
513       {
514        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
515       }
516       DataLazy_ptr lleft;
517       if (!left->isLazy())
518       {
519        lleft=DataLazy_ptr(new DataLazy(left));
520       }
521       else
522       {
523        lleft=dynamic_pointer_cast<DataLazy>(left);
524       }
525       m_readytype=lleft->m_readytype;
526       m_left=lleft;
527       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
528       m_samplesize=getNumDPPSample()*getNoValues();
529       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
530       m_children=m_left->m_children+1;
531       m_height=m_left->m_height+1;
532       SIZELIMIT
533    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
534    }
535    
536    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
537        : parent(left->getFunctionSpace(), left->getShape()),
538        m_op(op),
539        m_axis_offset(0),
540        m_transpose(0),
541        m_tol(tol)
542    {
543       if ((getOpgroup(op)!=G_UNARY_P))
544       {
545        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
546       }
547       DataLazy_ptr lleft;
548       if (!left->isLazy())
549       {
550        lleft=DataLazy_ptr(new DataLazy(left));
551       }
552       else
553       {
554        lleft=dynamic_pointer_cast<DataLazy>(left);
555       }
556       m_readytype=lleft->m_readytype;
557       m_left=lleft;
558       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
559       m_samplesize=getNumDPPSample()*getNoValues();
560       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
561       m_children=m_left->m_children+1;
562       m_height=m_left->m_height+1;
563       SIZELIMIT
564    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
565    }
566    
567  DataLazy::~DataLazy()  DataLazy::~DataLazy()
568  {  {
# Line 309  DataLazy::getBuffsRequired() const Line 576  DataLazy::getBuffsRequired() const
576  }  }
577    
578    
579    size_t
580    DataLazy::getMaxSampleSize() const
581    {
582        return m_maxsamplesize;
583    }
584    
585  /*  /*
586    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
587    This does the work for the collapse method.    This does the work for the collapse method.
# Line 328  DataLazy::collapseToReady() Line 601  DataLazy::collapseToReady()
601    DataReady_ptr pleft=m_left->collapseToReady();    DataReady_ptr pleft=m_left->collapseToReady();
602    Data left(pleft);    Data left(pleft);
603    Data right;    Data right;
604    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
605    {    {
606      right=Data(m_right->collapseToReady());      right=Data(m_right->collapseToReady());
607    }    }
# Line 427  DataLazy::collapseToReady() Line 700  DataLazy::collapseToReady()
700      case LEZ:      case LEZ:
701      result=left.whereNonPositive();      result=left.whereNonPositive();
702      break;      break;
703        case NEZ:
704        result=left.whereNonZero(m_tol);
705        break;
706        case EZ:
707        result=left.whereZero(m_tol);
708        break;
709        case SYM:
710        result=left.symmetric();
711        break;
712        case NSYM:
713        result=left.nonsymmetric();
714        break;
715        case PROD:
716        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
717        break;
718        case TRANS:
719        result=left.transpose(m_axis_offset);
720        break;
721        case TRACE:
722        result=left.trace(m_axis_offset);
723        break;
724      default:      default:
725      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)+".");
726    }    }
727    return result.borrowReadyPtr();    return result.borrowReadyPtr();
728  }  }
# Line 455  DataLazy::collapse() Line 749  DataLazy::collapse()
749  }  }
750    
751  /*  /*
752    \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.
753    \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.
754    \param v A vector to store intermediate results.    \param v A vector to store intermediate results.
755    \param offset Index in v to begin storing results.    \param offset Index in v to begin storing results.
# Line 483  DataLazy::resolveUnary(ValueType& v, siz Line 777  DataLazy::resolveUnary(ValueType& v, siz
777    switch (m_op)    switch (m_op)
778    {    {
779      case SIN:        case SIN:  
780      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
781      break;      break;
782      case COS:      case COS:
783      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
784      break;      break;
785      case TAN:      case TAN:
786      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
787      break;      break;
788      case ASIN:      case ASIN:
789      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
790      break;      break;
791      case ACOS:      case ACOS:
792      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
793      break;      break;
794      case ATAN:      case ATAN:
795      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
796      break;      break;
797      case SINH:      case SINH:
798      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
799      break;      break;
800      case COSH:      case COSH:
801      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
802      break;      break;
803      case TANH:      case TANH:
804      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
805      break;      break;
806      case ERF:      case ERF:
807  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
808      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
809  #else  #else
810      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
811      break;      break;
812  #endif  #endif
813     case ASINH:     case ASINH:
814  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
815      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
816  #else  #else
817      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
818  #endif    #endif  
819      break;      break;
820     case ACOSH:     case ACOSH:
821  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
822      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
823  #else  #else
824      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
825  #endif    #endif  
826      break;      break;
827     case ATANH:     case ATANH:
828  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
829      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
830  #else  #else
831      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
832  #endif    #endif  
833      break;      break;
834      case LOG10:      case LOG10:
835      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
836      break;      break;
837      case LOG:      case LOG:
838      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
839      break;      break;
840      case SIGN:      case SIGN:
841      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
842      break;      break;
843      case ABS:      case ABS:
844      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
845      break;      break;
846      case NEG:      case NEG:
847      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 558  DataLazy::resolveUnary(ValueType& v, siz Line 852  DataLazy::resolveUnary(ValueType& v, siz
852      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
853      break;      break;
854      case EXP:      case EXP:
855      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
856      break;      break;
857      case SQRT:      case SQRT:
858      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
859      break;      break;
860      case RECIP:      case RECIP:
861      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 578  DataLazy::resolveUnary(ValueType& v, siz Line 872  DataLazy::resolveUnary(ValueType& v, siz
872      case LEZ:      case LEZ:
873      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));
874      break;      break;
875    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
876        case NEZ:
877        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
878        break;
879        case EZ:
880        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
881        break;
882    
883      default:      default:
884      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 589  DataLazy::resolveUnary(ValueType& v, siz Line 890  DataLazy::resolveUnary(ValueType& v, siz
890    
891    
892    
893  #define PROC_OP(X) \  
894      for (int i=0;i<steps;++i,resultp+=resultStep) \  /*
895      { \    \brief Compute the value of the expression (unary operation) for the given sample.
896         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    \return Vector which stores the value of the subexpression for the given sample.
897         lroffset+=leftStep; \    \param v A vector to store intermediate results.
898         rroffset+=rightStep; \    \param offset Index in v to begin storing results.
899      \param sampleNo Sample number to evaluate.
900      \param roffset (output parameter) the offset in the return vector where the result begins.
901    
902      The return value will be an existing vector so do not deallocate it.
903      If the result is stored in v it should be stored at the offset given.
904      Everything from offset to the end of v should be considered available for this method to use.
905    */
906    DataTypes::ValueType*
907    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
908    {
909        // we assume that any collapsing has been done before we get here
910        // since we only have one argument we don't need to think about only
911        // processing single points.
912      if (m_readytype!='E')
913      {
914        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
915      }
916        // since we can't write the result over the input, we need a result offset further along
917      size_t subroffset=roffset+m_samplesize;
918    LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
919      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
920      roffset=offset;
921      size_t loop=0;
922      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
923      size_t step=getNoValues();
924      switch (m_op)
925      {
926        case SYM:
927        for (loop=0;loop<numsteps;++loop)
928        {
929            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
930            subroffset+=step;
931            offset+=step;
932        }
933        break;
934        case NSYM:
935        for (loop=0;loop<numsteps;++loop)
936        {
937            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
938            subroffset+=step;
939            offset+=step;
940        }
941        break;
942        default:
943        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
944      }
945      return &v;
946    }
947    
948    /*
949      \brief Compute the value of the expression (unary operation) for the given sample.
950      \return Vector which stores the value of the subexpression for the given sample.
951      \param v A vector to store intermediate results.
952      \param offset Index in v to begin storing results.
953      \param sampleNo Sample number to evaluate.
954      \param roffset (output parameter) the offset in the return vector where the result begins.
955    
956      The return value will be an existing vector so do not deallocate it.
957      If the result is stored in v it should be stored at the offset given.
958      Everything from offset to the end of v should be considered available for this method to use.
959    */
960    DataTypes::ValueType*
961    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
962    {
963        // we assume that any collapsing has been done before we get here
964        // since we only have one argument we don't need to think about only
965        // processing single points.
966      if (m_readytype!='E')
967      {
968        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
969      }
970        // since we can't write the result over the input, we need a result offset further along
971      size_t subroffset;
972      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
973    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
974    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
975      roffset=offset;
976      size_t loop=0;
977      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
978      size_t outstep=getNoValues();
979      size_t instep=m_left->getNoValues();
980    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
981      switch (m_op)
982      {
983        case TRACE:
984        for (loop=0;loop<numsteps;++loop)
985        {
986    size_t zz=sampleNo*getNumDPPSample()+loop;
987    if (zz==5800)
988    {
989    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
990    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
991    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
992    LAZYDEBUG(cerr << subroffset << endl;)
993    LAZYDEBUG(cerr << "output=" << offset << endl;)
994    }
995                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
996    if (zz==5800)
997    {
998    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
999    }
1000            subroffset+=instep;
1001            offset+=outstep;
1002        }
1003        break;
1004        case TRANS:
1005        for (loop=0;loop<numsteps;++loop)
1006        {
1007                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1008            subroffset+=instep;
1009            offset+=outstep;
1010        }
1011        break;
1012        default:
1013        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1014      }
1015      return &v;
1016    }
1017    
1018    
1019    #define PROC_OP(TYPE,X)                               \
1020        for (int j=0;j<onumsteps;++j)\
1021        {\
1022          for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1023          { \
1024    LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1025    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1026             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1027    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1028             lroffset+=leftstep; \
1029             rroffset+=rightstep; \
1030          }\
1031          lroffset+=oleftstep;\
1032          rroffset+=orightstep;\
1033      }      }
1034    
1035  /*  /*
# Line 621  DataLazy::resolveUnary(ValueType& v, siz Line 1056  DataLazy::resolveUnary(ValueType& v, siz
1056  DataTypes::ValueType*  DataTypes::ValueType*
1057  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
1058  {  {
1059  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1060    
1061    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
1062      // first work out which of the children are expanded      // first work out which of the children are expanded
1063    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1064    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1065    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1066    int steps=(bigloops?1:getNumDPPSample());    {
1067    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'.");
1068    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1069    {    bool leftScalar=(m_left->getRank()==0);
1070      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1071      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1072      chunksize=1;    // for scalar    {
1073    }          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1074    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1075    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1076    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1077      size_t chunksize=1;           // how many doubles will be processed in one go
1078      int leftstep=0;       // how far should the left offset advance after each step
1079      int rightstep=0;
1080      int numsteps=0;       // total number of steps for the inner loop
1081      int oleftstep=0;  // the o variables refer to the outer loop
1082      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1083      int onumsteps=1;
1084      
1085      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1086      bool RES=(rightExp && rightScalar);
1087      bool LS=(!leftExp && leftScalar); // left is a single scalar
1088      bool RS=(!rightExp && rightScalar);
1089      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1090      bool RN=(!rightExp && !rightScalar);
1091      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1092      bool REN=(rightExp && !rightScalar);
1093    
1094      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1095      {
1096        chunksize=m_left->getNumDPPSample()*leftsize;
1097        leftstep=0;
1098        rightstep=0;
1099        numsteps=1;
1100      }
1101      else if (LES || RES)
1102      {
1103        chunksize=1;
1104        if (LES)        // left is an expanded scalar
1105        {
1106            if (RS)
1107            {
1108               leftstep=1;
1109               rightstep=0;
1110               numsteps=m_left->getNumDPPSample();
1111            }
1112            else        // RN or REN
1113            {
1114               leftstep=0;
1115               oleftstep=1;
1116               rightstep=1;
1117               orightstep=(RN ? -(int)rightsize : 0);
1118               numsteps=rightsize;
1119               onumsteps=m_left->getNumDPPSample();
1120            }
1121        }
1122        else        // right is an expanded scalar
1123        {
1124            if (LS)
1125            {
1126               rightstep=1;
1127               leftstep=0;
1128               numsteps=m_right->getNumDPPSample();
1129            }
1130            else
1131            {
1132               rightstep=0;
1133               orightstep=1;
1134               leftstep=1;
1135               oleftstep=(LN ? -(int)leftsize : 0);
1136               numsteps=leftsize;
1137               onumsteps=m_right->getNumDPPSample();
1138            }
1139        }
1140      }
1141      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1142      {
1143        if (LEN)    // and Right will be a single value
1144        {
1145            chunksize=rightsize;
1146            leftstep=rightsize;
1147            rightstep=0;
1148            numsteps=m_left->getNumDPPSample();
1149            if (RS)
1150            {
1151               numsteps*=leftsize;
1152            }
1153        }
1154        else    // REN
1155        {
1156            chunksize=leftsize;
1157            rightstep=leftsize;
1158            leftstep=0;
1159            numsteps=m_right->getNumDPPSample();
1160            if (LS)
1161            {
1162               numsteps*=rightsize;
1163            }
1164        }
1165      }
1166    
1167      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1168      // Get the values of sub-expressions      // Get the values of sub-expressions
1169    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1170    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
1171      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1172      // the right child starts further along.      // the right child starts further along.
1173    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1174    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1175    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1176    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1177    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1178    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1179    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1180    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
1181    switch(m_op)    switch(m_op)
1182    {    {
1183      case ADD:      case ADD:
1184      PROC_OP(plus<double>());          PROC_OP(NO_ARG,plus<double>());
1185      break;      break;
1186      case SUB:      case SUB:
1187      PROC_OP(minus<double>());      PROC_OP(NO_ARG,minus<double>());
1188      break;      break;
1189      case MUL:      case MUL:
1190      PROC_OP(multiplies<double>());      PROC_OP(NO_ARG,multiplies<double>());
1191      break;      break;
1192      case DIV:      case DIV:
1193      PROC_OP(divides<double>());      PROC_OP(NO_ARG,divides<double>());
1194      break;      break;
1195      case POW:      case POW:
1196      PROC_OP(::pow);         PROC_OP(double (double,double),::pow);
1197      break;      break;
1198      default:      default:
1199      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
# Line 671  cout << "Resolve binary: " << toString() Line 1205  cout << "Resolve binary: " << toString()
1205    
1206    
1207  /*  /*
1208      \brief Compute the value of the expression (tensor product) for the given sample.
1209      \return Vector which stores the value of the subexpression for the given sample.
1210      \param v A vector to store intermediate results.
1211      \param offset Index in v to begin storing results.
1212      \param sampleNo Sample number to evaluate.
1213      \param roffset (output parameter) the offset in the return vector where the result begins.
1214    
1215      The return value will be an existing vector so do not deallocate it.
1216      If the result is stored in v it should be stored at the offset given.
1217      Everything from offset to the end of v should be considered available for this method to use.
1218    */
1219    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1220    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1221    // unlike the other resolve helpers, we must treat these datapoints separately.
1222    DataTypes::ValueType*
1223    DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1224    {
1225    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1226    
1227      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1228        // first work out which of the children are expanded
1229      bool leftExp=(m_left->m_readytype=='E');
1230      bool rightExp=(m_right->m_readytype=='E');
1231      int steps=getNumDPPSample();
1232      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1233      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
1234      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
1235        // Get the values of sub-expressions (leave a gap of one sample for the result).
1236      int gap=offset+m_left->getMaxSampleSize();    // actually only needs to be m_left->m_samplesize
1237      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1238      gap+=m_right->getMaxSampleSize();
1239      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1240    LAZYDEBUG(cout << "Post sub calls: " << toString() << endl;)
1241    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1242    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1243    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1244    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1245    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1246      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1247      switch(m_op)
1248      {
1249        case PROD:
1250        for (int i=0;i<steps;++i,resultp+=resultStep)
1251        {
1252    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1253    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1254    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1255              const double *ptr_0 = &((*left)[lroffset]);
1256              const double *ptr_1 = &((*right)[rroffset]);
1257    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1258    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1259              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1260          lroffset+=leftStep;
1261          rroffset+=rightStep;
1262        }
1263        break;
1264        default:
1265        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1266      }
1267      roffset=offset;
1268      return &v;
1269    }
1270    
1271    
1272    
1273    /*
1274    \brief Compute the value of the expression for the given sample.    \brief Compute the value of the expression for the given sample.
1275    \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.
1276    \param v A vector to store intermediate results.    \param v A vector to store intermediate results.
# Line 689  cout << "Resolve binary: " << toString() Line 1289  cout << "Resolve binary: " << toString()
1289  const DataTypes::ValueType*  const DataTypes::ValueType*
1290  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1291  {  {
1292  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1293      // 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
1294    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1295    {    {
# Line 701  cout << "Resolve sample " << toString() Line 1301  cout << "Resolve sample " << toString()
1301      if (m_readytype=='C')      if (m_readytype=='C')
1302      {      {
1303      roffset=0;      roffset=0;
1304    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1305      return &(vec);      return &(vec);
1306      }      }
1307      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1308    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1309      return &(vec);      return &(vec);
1310    }    }
1311    if (m_readytype!='E')    if (m_readytype!='E')
# Line 712  cout << "Resolve sample " << toString() Line 1314  cout << "Resolve sample " << toString()
1314    }    }
1315    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1316    {    {
1317    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY:
1318      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1319    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1320      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1321      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1322      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1323    default:    default:
1324      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)+".");
1325    }    }
1326    
1327    }
1328    
1329    // This needs to do the work of the idenity constructor
1330    void
1331    DataLazy::resolveToIdentity()
1332    {
1333       if (m_op==IDENTITY)
1334        return;
1335       DataReady_ptr p=resolve();
1336       makeIdentity(p);
1337  }  }
1338    
1339    void DataLazy::makeIdentity(const DataReady_ptr& p)
1340    {
1341       m_op=IDENTITY;
1342       m_axis_offset=0;
1343       m_transpose=0;
1344       m_SL=m_SM=m_SR=0;
1345       m_children=m_height=0;
1346       m_id=p;
1347       if(p->isConstant()) {m_readytype='C';}
1348       else if(p->isExpanded()) {m_readytype='E';}
1349       else if (p->isTagged()) {m_readytype='T';}
1350       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1351       m_buffsRequired=1;
1352       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1353       m_maxsamplesize=m_samplesize;
1354       m_left.reset();
1355       m_right.reset();
1356    }
1357    
1358  // 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.
1359  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
# Line 726  DataReady_ptr Line 1361  DataReady_ptr
1361  DataLazy::resolve()  DataLazy::resolve()
1362  {  {
1363    
1364  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1365  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1366    
1367    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
1368    {    {
# Line 738  cout << "Buffers=" << m_buffsRequired << Line 1373  cout << "Buffers=" << m_buffsRequired <<
1373      return m_id;      return m_id;
1374    }    }
1375      // 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'
1376    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
1377      // storage to evaluate its expression      // storage to evaluate its expression
1378    int numthreads=1;    int numthreads=1;
1379  #ifdef _OPENMP  #ifdef _OPENMP
1380    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
   int threadnum=0;  
1381  #endif  #endif
1382    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1383  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1384    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1385    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1386    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
# Line 755  cout << "Buffer created with size=" << v Line 1389  cout << "Buffer created with size=" << v
1389    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1390    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1391    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1392    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1393      #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1394    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1395    {    {
1396  cout << "################################# " << sample << endl;        if (sample==0)  {ENABLEDEBUG}
1397    
1398    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1399    LAZYDEBUG(cout << "################################# " << sample << endl;)
1400  #ifdef _OPENMP  #ifdef _OPENMP
1401      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1402  #else  #else
1403      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.
1404  #endif  #endif
1405  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1406    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1407      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1408  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1409      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
1410      {      {
1411    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1412      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1413      }      }
1414  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1415    LAZYDEBUG(cerr << "*********************************" << endl;)
1416        DISABLEDEBUG
1417    }    }
1418    return resptr;    return resptr;
1419  }  }
# Line 789  DataLazy::toString() const Line 1431  DataLazy::toString() const
1431  void  void
1432  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1433  {  {
1434      oss << "[" << m_children <<";"<<m_height <<"]";
1435    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1436    {    {
1437    case G_IDENTITY:    case G_IDENTITY:
# Line 818  DataLazy::intoString(ostringstream& oss) Line 1461  DataLazy::intoString(ostringstream& oss)
1461      oss << ')';      oss << ')';
1462      break;      break;
1463    case G_UNARY:    case G_UNARY:
1464      case G_UNARY_P:
1465      case G_NP1OUT:
1466      case G_NP1OUT_P:
1467      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1468      m_left->intoString(oss);      m_left->intoString(oss);
1469      oss << ')';      oss << ')';
1470      break;      break;
1471      case G_TENSORPROD:
1472        oss << opToString(m_op) << '(';
1473        m_left->intoString(oss);
1474        oss << ", ";
1475        m_right->intoString(oss);
1476        oss << ')';
1477        break;
1478    default:    default:
1479      oss << "UNKNOWN";      oss << "UNKNOWN";
1480    }    }
1481  }  }
1482    
 // 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.  
1483  DataAbstract*  DataAbstract*
1484  DataLazy::deepCopy()  DataLazy::deepCopy()
1485  {  {
1486    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
1487    {    {
1488      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1489      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1490      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1491      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1492      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1493      default:
1494        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1495    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
1496  }  }
1497    
1498    
1499    // There is no single, natural interpretation of getLength on DataLazy.
1500    // Instances of DataReady can look at the size of their vectors.
1501    // For lazy though, it could be the size the data would be if it were resolved;
1502    // or it could be some function of the lengths of the DataReady instances which
1503    // form part of the expression.
1504    // Rather than have people making assumptions, I have disabled the method.
1505  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1506  DataLazy::getLength() const  DataLazy::getLength() const
1507  {  {
1508    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1509  }  }
1510    
1511    
# Line 853  DataLazy::getSlice(const DataTypes::Regi Line 1515  DataLazy::getSlice(const DataTypes::Regi
1515    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
1516  }  }
1517    
1518    
1519    // To do this we need to rely on our child nodes
1520    DataTypes::ValueType::size_type
1521    DataLazy::getPointOffset(int sampleNo,
1522                     int dataPointNo)
1523    {
1524      if (m_op==IDENTITY)
1525      {
1526        return m_id->getPointOffset(sampleNo,dataPointNo);
1527      }
1528      if (m_readytype!='E')
1529      {
1530        collapse();
1531        return m_id->getPointOffset(sampleNo,dataPointNo);
1532      }
1533      // at this point we do not have an identity node and the expression will be Expanded
1534      // so we only need to know which child to ask
1535      if (m_left->m_readytype=='E')
1536      {
1537        return m_left->getPointOffset(sampleNo,dataPointNo);
1538      }
1539      else
1540      {
1541        return m_right->getPointOffset(sampleNo,dataPointNo);
1542      }
1543    }
1544    
1545    // To do this we need to rely on our child nodes
1546  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1547  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1548                   int dataPointNo) const                   int dataPointNo) const
1549  {  {
1550    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
1551      {
1552        return m_id->getPointOffset(sampleNo,dataPointNo);
1553      }
1554      if (m_readytype=='E')
1555      {
1556        // at this point we do not have an identity node and the expression will be Expanded
1557        // so we only need to know which child to ask
1558        if (m_left->m_readytype=='E')
1559        {
1560        return m_left->getPointOffset(sampleNo,dataPointNo);
1561        }
1562        else
1563        {
1564        return m_right->getPointOffset(sampleNo,dataPointNo);
1565        }
1566      }
1567      if (m_readytype=='C')
1568      {
1569        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1570      }
1571      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1572  }  }
1573    
1574  // It would seem that DataTagged will need to be treated differently since even after setting all tags  // It would seem that DataTagged will need to be treated differently since even after setting all tags

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

  ViewVC Help
Powered by ViewVC 1.1.26