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

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

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

branches/schroedinger_upto1946/escript/src/DataLazy.cpp revision 1972 by jfenwick, Thu Nov 6 01:49:39 2008 UTC trunk/escript/src/DataLazy.cpp revision 2082 by caltinay, Fri Nov 21 01:46:05 2008 UTC
# Line 48  I will refer to individual DataLazy obje Line 48  I will refer to individual DataLazy obje
48  Each node also stores:  Each node also stores:
49  - 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
50      evaluated.      evaluated.
 - m_length ~ how many values would be stored in the answer if the expression was evaluated.  
51  - 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
52      evaluate the expression.      evaluate the expression.
53  - m_samplesize ~ the number of doubles stored in a sample.  - m_samplesize ~ the number of doubles stored in a sample.
# Line 70  The convention that I use, is that the r Line 69  The convention that I use, is that the r
69    
70  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.
71  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.
72    
73    To add a new operator you need to do the following (plus anything I might have forgotten):
74    1) Add to the ES_optype.
75    2) determine what opgroup your operation belongs to (X)
76    3) add a string for the op to the end of ES_opstrings
77    4) increase ES_opcount
78    5) add an entry (X) to opgroups
79    6) add an entry to the switch in collapseToReady
80    7) add an entry to resolveX
81  */  */
82    
83    
# Line 79  using namespace boost; Line 87  using namespace boost;
87  namespace escript  namespace escript
88  {  {
89    
 const std::string&  
 opToString(ES_optype op);  
   
90  namespace  namespace
91  {  {
92    
# Line 90  enum ES_opgroup Line 95  enum ES_opgroup
95     G_UNKNOWN,     G_UNKNOWN,
96     G_IDENTITY,     G_IDENTITY,
97     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
98     G_UNARY      // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
99       G_NP1OUT,        // non-pointwise op with one output
100       G_TENSORPROD     // general tensor product
101  };  };
102    
103    
# Line 101  string ES_opstrings[]={"UNKNOWN","IDENTI Line 108  string ES_opstrings[]={"UNKNOWN","IDENTI
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=33;              "symmetric","nonsymmetric",
113                "prod"};
114    int ES_opcount=36;
115  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,
116              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
117              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
118              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
119              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
120              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,            // 33
121                G_NP1OUT,G_NP1OUT,
122                G_TENSORPROD};
123  inline  inline
124  ES_opgroup  ES_opgroup
125  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 143  resultFS(DataAbstract_ptr left, DataAbst Line 154  resultFS(DataAbstract_ptr left, DataAbst
154  }  }
155    
156  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
157    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
158  DataTypes::ShapeType  DataTypes::ShapeType
159  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
160  {  {
161      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
162      {      {
163        if (getOpgroup(op)!=G_BINARY)        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
164        {        {
165          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.");
166        }        }
# Line 165  resultShape(DataAbstract_ptr left, DataA Line 177  resultShape(DataAbstract_ptr left, DataA
177      return left->getShape();      return left->getShape();
178  }  }
179    
180  // determine the number of points in the result of "left op right"  // determine the output shape for the general tensor product operation
181  size_t  // the additional parameters return information required later for the product
182  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  // the majority of this code is copy pasted from C_General_Tensor_Product
183    DataTypes::ShapeType
184    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
185  {  {
186     switch (getOpgroup(op))      
187     {    // Get rank and shape of inputs
188     case G_BINARY: return left->getLength();    int rank0 = left->getRank();
189     case G_UNARY: return left->getLength();    int rank1 = right->getRank();
190     default:    const DataTypes::ShapeType& shape0 = left->getShape();
191      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");    const DataTypes::ShapeType& shape1 = right->getShape();
192     }  
193      // Prepare for the loops of the product and verify compatibility of shapes
194      int start0=0, start1=0;
195      if (transpose == 0)       {}
196      else if (transpose == 1)  { start0 = axis_offset; }
197      else if (transpose == 2)  { start1 = rank1-axis_offset; }
198      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
199    
200    
201      // Adjust the shapes for transpose
202      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
203      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
204      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
205      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
206    
207      // Prepare for the loops of the product
208      SL=1, SM=1, SR=1;
209      for (int i=0; i<rank0-axis_offset; i++)   {
210        SL *= tmpShape0[i];
211      }
212      for (int i=rank0-axis_offset; i<rank0; i++)   {
213        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
214          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
215        }
216        SM *= tmpShape0[i];
217      }
218      for (int i=axis_offset; i<rank1; i++)     {
219        SR *= tmpShape1[i];
220      }
221    
222      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
223      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
224      {         // block to limit the scope of out_index
225         int out_index=0;
226         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
227         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
228      }
229      return shape2;
230  }  }
231    
232    
233    // determine the number of points in the result of "left op right"
234    // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here
235    // size_t
236    // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
237    // {
238    //    switch (getOpgroup(op))
239    //    {
240    //    case G_BINARY: return left->getLength();
241    //    case G_UNARY: return left->getLength();
242    //    case G_NP1OUT: return left->getLength();
243    //    default:
244    //  throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
245    //    }
246    // }
247    
248  // 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
249    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
250    // The same goes for G_TENSORPROD
251  int  int
252  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
253  {  {
# Line 187  calcBuffs(const DataLazy_ptr& left, cons Line 256  calcBuffs(const DataLazy_ptr& left, cons
256     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
257     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
258     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY: return max(left->getBuffsRequired(),1);
259       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
260       case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
261     default:     default:
262      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
263     }     }
# Line 211  opToString(ES_optype op) Line 282  opToString(ES_optype op)
282    
283  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
284      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape()),
285      m_op(IDENTITY)      m_op(IDENTITY),
286        m_axis_offset(0),
287        m_transpose(0),
288        m_SL(0), m_SM(0), m_SR(0)
289  {  {
290     if (p->isLazy())     if (p->isLazy())
291     {     {
# Line 228  DataLazy::DataLazy(DataAbstract_ptr p) Line 302  DataLazy::DataLazy(DataAbstract_ptr p)
302      else if (p->isTagged()) {m_readytype='T';}      else if (p->isTagged()) {m_readytype='T';}
303      else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}      else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
304     }     }
    m_length=p->getLength();  
305     m_buffsRequired=1;     m_buffsRequired=1;
306     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
307       m_maxsamplesize=m_samplesize;
308  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
309  }  }
310    
# Line 239  cout << "(1)Lazy created with " << m_sam Line 313  cout << "(1)Lazy created with " << m_sam
313    
314  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
315      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
316      m_op(op)      m_op(op),
317        m_axis_offset(0),
318        m_transpose(0),
319        m_SL(0), m_SM(0), m_SR(0)
320  {  {
321     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
322     {     {
323      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.");
324     }     }
325    
326     DataLazy_ptr lleft;     DataLazy_ptr lleft;
327     if (!left->isLazy())     if (!left->isLazy())
328     {     {
# Line 255  DataLazy::DataLazy(DataAbstract_ptr left Line 333  DataLazy::DataLazy(DataAbstract_ptr left
333      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
334     }     }
335     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
    m_length=left->getLength();  
336     m_left=lleft;     m_left=lleft;
337     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
338     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
339       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
340  }  }
341    
342    
343  // In this constructor we need to consider interpolation  // In this constructor we need to consider interpolation
344  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
345      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
346      m_op(op)      m_op(op),
347        m_SL(0), m_SM(0), m_SR(0)
348  {  {
349     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_BINARY))
350     {     {
351      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.");
352     }     }
# Line 279  DataLazy::DataLazy(DataAbstract_ptr left Line 358  DataLazy::DataLazy(DataAbstract_ptr left
358      Data tmp(ltemp,fs);      Data tmp(ltemp,fs);
359      left=tmp.borrowDataPtr();      left=tmp.borrowDataPtr();
360     }     }
361     if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
362     {     {
363      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
364      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
# Line 316  DataLazy::DataLazy(DataAbstract_ptr left Line 395  DataLazy::DataLazy(DataAbstract_ptr left
395     {     {
396      m_readytype='C';      m_readytype='C';
397     }     }
398     m_length=resultLength(m_left,m_right,m_op);     m_samplesize=getNumDPPSample()*getNoValues();
399     m_samplesize=getNumDPPSample()*getNoValues();         m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
400     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
401  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
402  }  }
403    
404    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
405        : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
406        m_op(op),
407        m_axis_offset(axis_offset),
408        m_transpose(transpose)
409    {
410       if ((getOpgroup(op)!=G_TENSORPROD))
411       {
412        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
413       }
414       if ((transpose>2) || (transpose<0))
415       {
416        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
417       }
418       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
419       {
420        FunctionSpace fs=getFunctionSpace();
421        Data ltemp(left);
422        Data tmp(ltemp,fs);
423        left=tmp.borrowDataPtr();
424       }
425       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
426       {
427        Data tmp(Data(right),getFunctionSpace());
428        right=tmp.borrowDataPtr();
429       }
430       left->operandCheck(*right);
431    
432       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
433       {
434        m_left=dynamic_pointer_cast<DataLazy>(left);
435       }
436       else
437       {
438        m_left=DataLazy_ptr(new DataLazy(left));
439       }
440       if (right->isLazy())
441       {
442        m_right=dynamic_pointer_cast<DataLazy>(right);
443       }
444       else
445       {
446        m_right=DataLazy_ptr(new DataLazy(right));
447       }
448       char lt=m_left->m_readytype;
449       char rt=m_right->m_readytype;
450       if (lt=='E' || rt=='E')
451       {
452        m_readytype='E';
453       }
454       else if (lt=='T' || rt=='T')
455       {
456        m_readytype='T';
457       }
458       else
459       {
460        m_readytype='C';
461       }
462       m_samplesize=getNumDPPSample()*getNoValues();
463       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
464       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
465    cout << "(4)Lazy created with " << m_samplesize << endl;
466    }
467    
468    
469  DataLazy::~DataLazy()  DataLazy::~DataLazy()
470  {  {
# Line 335  DataLazy::getBuffsRequired() const Line 478  DataLazy::getBuffsRequired() const
478  }  }
479    
480    
481    size_t
482    DataLazy::getMaxSampleSize() const
483    {
484        return m_maxsamplesize;
485    }
486    
487  /*  /*
488    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
489    This does the work for the collapse method.    This does the work for the collapse method.
# Line 354  DataLazy::collapseToReady() Line 503  DataLazy::collapseToReady()
503    DataReady_ptr pleft=m_left->collapseToReady();    DataReady_ptr pleft=m_left->collapseToReady();
504    Data left(pleft);    Data left(pleft);
505    Data right;    Data right;
506    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
507    {    {
508      right=Data(m_right->collapseToReady());      right=Data(m_right->collapseToReady());
509    }    }
# Line 453  DataLazy::collapseToReady() Line 602  DataLazy::collapseToReady()
602      case LEZ:      case LEZ:
603      result=left.whereNonPositive();      result=left.whereNonPositive();
604      break;      break;
605        case SYM:
606        result=left.symmetric();
607        break;
608        case NSYM:
609        result=left.nonsymmetric();
610        break;
611        case PROD:
612        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
613        break;
614      default:      default:
615      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)+".");
616    }    }
617    return result.borrowReadyPtr();    return result.borrowReadyPtr();
618  }  }
# Line 481  DataLazy::collapse() Line 639  DataLazy::collapse()
639  }  }
640    
641  /*  /*
642    \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.
643    \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.
644    \param v A vector to store intermediate results.    \param v A vector to store intermediate results.
645    \param offset Index in v to begin storing results.    \param offset Index in v to begin storing results.
# Line 536  DataLazy::resolveUnary(ValueType& v, siz Line 694  DataLazy::resolveUnary(ValueType& v, siz
694      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
695      break;      break;
696      case ERF:      case ERF:
697  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
698      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
699  #else  #else
700      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
701      break;      break;
702  #endif  #endif
703     case ASINH:     case ASINH:
704  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
705      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
706  #else  #else
707      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
708  #endif    #endif  
709      break;      break;
710     case ACOSH:     case ACOSH:
711  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
712      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
713  #else  #else
714      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
715  #endif    #endif  
716      break;      break;
717     case ATANH:     case ATANH:
718  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
719      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
720  #else  #else
721      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
# Line 612  DataLazy::resolveUnary(ValueType& v, siz Line 770  DataLazy::resolveUnary(ValueType& v, siz
770  }  }
771    
772    
773    /*
774      \brief Compute the value of the expression (unary operation) for the given sample.
775      \return Vector which stores the value of the subexpression for the given sample.
776      \param v A vector to store intermediate results.
777      \param offset Index in v to begin storing results.
778      \param sampleNo Sample number to evaluate.
779      \param roffset (output parameter) the offset in the return vector where the result begins.
780    
781      The return value will be an existing vector so do not deallocate it.
782      If the result is stored in v it should be stored at the offset given.
783      Everything from offset to the end of v should be considered available for this method to use.
784    */
785    DataTypes::ValueType*
786    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
787    {
788        // we assume that any collapsing has been done before we get here
789        // since we only have one argument we don't need to think about only
790        // processing single points.
791      if (m_readytype!='E')
792      {
793        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
794      }
795        // since we can't write the result over the input, we need a result offset further along
796      size_t subroffset=roffset+m_samplesize;
797      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
798      roffset=offset;
799      switch (m_op)
800      {
801        case SYM:
802        DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
803        break;
804        case NSYM:
805        DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
806        break;
807        default:
808        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
809      }
810      return &v;
811    }
812    
813    
814    
815    
816  #define PROC_OP(X) \  #define PROC_OP(TYPE,X)                               \
817      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int i=0;i<steps;++i,resultp+=resultStep) \
818      { \      { \
819         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
820         lroffset+=leftStep; \         lroffset+=leftStep; \
821         rroffset+=rightStep; \         rroffset+=rightStep; \
822      }      }
# Line 673  cout << "Resolve binary: " << toString() Line 871  cout << "Resolve binary: " << toString()
871    switch(m_op)    switch(m_op)
872    {    {
873      case ADD:      case ADD:
874      PROC_OP(plus<double>());          PROC_OP(NO_ARG,plus<double>());
875      break;      break;
876      case SUB:      case SUB:
877      PROC_OP(minus<double>());      PROC_OP(NO_ARG,minus<double>());
878      break;      break;
879      case MUL:      case MUL:
880      PROC_OP(multiplies<double>());      PROC_OP(NO_ARG,multiplies<double>());
881      break;      break;
882      case DIV:      case DIV:
883      PROC_OP(divides<double>());      PROC_OP(NO_ARG,divides<double>());
884      break;      break;
885      case POW:      case POW:
886      PROC_OP(::pow);         PROC_OP(double (double,double),::pow);
887      break;      break;
888      default:      default:
889      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 695  cout << "Resolve binary: " << toString() Line 893  cout << "Resolve binary: " << toString()
893  }  }
894    
895    
896    /*
897      \brief Compute the value of the expression (tensor product) for the given sample.
898      \return Vector which stores the value of the subexpression for the given sample.
899      \param v A vector to store intermediate results.
900      \param offset Index in v to begin storing results.
901      \param sampleNo Sample number to evaluate.
902      \param roffset (output parameter) the offset in the return vector where the result begins.
903    
904      The return value will be an existing vector so do not deallocate it.
905      If the result is stored in v it should be stored at the offset given.
906      Everything from offset to the end of v should be considered available for this method to use.
907    */
908    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
909    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
910    // unlike the other resolve helpers, we must treat these datapoints separately.
911    DataTypes::ValueType*
912    DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
913    {
914    cout << "Resolve TensorProduct: " << toString() << endl;
915    
916      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
917        // first work out which of the children are expanded
918      bool leftExp=(m_left->m_readytype=='E');
919      bool rightExp=(m_right->m_readytype=='E');
920      int steps=getNumDPPSample();
921      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
922      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
923      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
924        // Get the values of sub-expressions (leave a gap of one sample for the result).
925      const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);
926      const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);
927      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
928      switch(m_op)
929      {
930        case PROD:
931        for (int i=0;i<steps;++i,resultp+=resultStep)
932        {
933              const double *ptr_0 = &((*left)[lroffset]);
934              const double *ptr_1 = &((*right)[rroffset]);
935              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
936          lroffset+=leftStep;
937          rroffset+=rightStep;
938        }
939        break;
940        default:
941        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
942      }
943      roffset=offset;
944      return &v;
945    }
946    
947    
948    
949  /*  /*
950    \brief Compute the value of the expression for the given sample.    \brief Compute the value of the expression for the given sample.
# Line 740  cout << "Resolve sample " << toString() Line 990  cout << "Resolve sample " << toString()
990    {    {
991    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
992    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
993      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
994      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
995    default:    default:
996      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)+".");
997    }    }
# Line 764  cout << "Buffers=" << m_buffsRequired << Line 1016  cout << "Buffers=" << m_buffsRequired <<
1016      return m_id;      return m_id;
1017    }    }
1018      // 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'
1019    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
1020      // storage to evaluate its expression      // storage to evaluate its expression
1021    int numthreads=1;    int numthreads=1;
1022  #ifdef _OPENMP  #ifdef _OPENMP
1023    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
   int threadnum=0;  
1024  #endif  #endif
1025    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1026  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
# Line 781  cout << "Buffer created with size=" << v Line 1032  cout << "Buffer created with size=" << v
1032    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1033    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1034    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1035    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1036    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1037    {    {
1038  cout << "################################# " << sample << endl;  cout << "################################# " << sample << endl;
# Line 844  DataLazy::intoString(ostringstream& oss) Line 1095  DataLazy::intoString(ostringstream& oss)
1095      oss << ')';      oss << ')';
1096      break;      break;
1097    case G_UNARY:    case G_UNARY:
1098      case G_NP1OUT:
1099      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1100      m_left->intoString(oss);      m_left->intoString(oss);
1101      oss << ')';      oss << ')';
1102      break;      break;
1103      case G_TENSORPROD:
1104        oss << opToString(m_op) << '(';
1105        m_left->intoString(oss);
1106        oss << ", ";
1107        m_right->intoString(oss);
1108        oss << ')';
1109        break;
1110    default:    default:
1111      oss << "UNKNOWN";      oss << "UNKNOWN";
1112    }    }
# Line 861  DataLazy::deepCopy() Line 1120  DataLazy::deepCopy()
1120    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1121    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1122    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1123      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1124      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1125    default:    default:
1126      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1127    }    }
1128  }  }
1129    
1130    
1131    // There is no single, natural interpretation of getLength on DataLazy.
1132    // Instances of DataReady can look at the size of their vectors.
1133    // For lazy though, it could be the size the data would be if it were resolved;
1134    // or it could be some function of the lengths of the DataReady instances which
1135    // form part of the expression.
1136    // Rather than have people making assumptions, I have disabled the method.
1137  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1138  DataLazy::getLength() const  DataLazy::getLength() const
1139  {  {
1140    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1141  }  }
1142    
1143    

Legend:
Removed from v.1972  
changed lines
  Added in v.2082

  ViewVC Help
Powered by ViewVC 1.1.26