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

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

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1943 by jfenwick, Wed Oct 29 04:05:14 2008 UTC trunk/escript/src/DataLazy.cpp revision 2050 by phornby, Mon Nov 17 08:59:57 2008 UTC
# Line 26  Line 26 
26  #include "DataTypes.h"  #include "DataTypes.h"
27  #include "Data.h"  #include "Data.h"
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29    #include "Utils.h"
30    
31  /*  /*
32  How does DataLazy work?  How does DataLazy work?
# Line 69  The convention that I use, is that the r Line 70  The convention that I use, is that the r
70    
71  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.
72  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.
73    
74    To add a new operator you need to do the following (plus anything I might have forgotten):
75    1) Add to the ES_optype.
76    2) determine what opgroup your operation belongs to (X)
77    3) add a string for the op to the end of ES_opstrings
78    4) increase ES_opcount
79    5) add an entry (X) to opgroups
80    6) add an entry to the switch in collapseToReady
81    7) add an entry to resolveX
82  */  */
83    
84    
# Line 78  using namespace boost; Line 88  using namespace boost;
88  namespace escript  namespace escript
89  {  {
90    
 const std::string&  
 opToString(ES_optype op);  
   
91  namespace  namespace
92  {  {
93    
# Line 89  enum ES_opgroup Line 96  enum ES_opgroup
96     G_UNKNOWN,     G_UNKNOWN,
97     G_IDENTITY,     G_IDENTITY,
98     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
99     G_UNARY      // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
100       G_NP1OUT     // non-pointwise op with one output
101  };  };
102    
103    
# Line 100  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    int ES_opcount=35;
114  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,
115              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
116              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
117              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
118              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
119              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,            // 33
120                G_NP1OUT,G_NP1OUT};
121  inline  inline
122  ES_opgroup  ES_opgroup
123  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 147  resultShape(DataAbstract_ptr left, DataA Line 157  resultShape(DataAbstract_ptr left, DataA
157  {  {
158      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
159      {      {
160        if (getOpgroup(op)!=G_BINARY)        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
161        {        {
162          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.");
163        }        }
# Line 172  resultLength(DataAbstract_ptr left, Data Line 182  resultLength(DataAbstract_ptr left, Data
182     {     {
183     case G_BINARY: return left->getLength();     case G_BINARY: return left->getLength();
184     case G_UNARY: return left->getLength();     case G_UNARY: return left->getLength();
185       case G_NP1OUT: return left->getLength();
186     default:     default:
187      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
188     }     }
189  }  }
190    
191  // determine the number of samples requires to evaluate an expression combining left and right  // determine the number of samples requires to evaluate an expression combining left and right
192    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
193  int  int
194  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
195  {  {
# Line 186  calcBuffs(const DataLazy_ptr& left, cons Line 198  calcBuffs(const DataLazy_ptr& left, cons
198     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
199     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
200     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY: return max(left->getBuffsRequired(),1);
201       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
202     default:     default:
203      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
204     }     }
# Line 240  DataLazy::DataLazy(DataAbstract_ptr left Line 253  DataLazy::DataLazy(DataAbstract_ptr left
253      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
254      m_op(op)      m_op(op)
255  {  {
256     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
257     {     {
258      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
259     }     }
# Line 256  DataLazy::DataLazy(DataAbstract_ptr left Line 269  DataLazy::DataLazy(DataAbstract_ptr left
269     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
270     m_length=left->getLength();     m_length=left->getLength();
271     m_left=lleft;     m_left=lleft;
272     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
273     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
274  }  }
275    
# Line 266  DataLazy::DataLazy(DataAbstract_ptr left Line 279  DataLazy::DataLazy(DataAbstract_ptr left
279      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
280      m_op(op)      m_op(op)
281  {  {
282     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_BINARY))
283     {     {
284    cout << opToString(op) << endl;
285      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
286     }     }
287    
# Line 452  DataLazy::collapseToReady() Line 466  DataLazy::collapseToReady()
466      case LEZ:      case LEZ:
467      result=left.whereNonPositive();      result=left.whereNonPositive();
468      break;      break;
469        case SYM:
470        result=left.symmetric();
471        break;
472        case NSYM:
473        result=left.nonsymmetric();
474        break;
475      default:      default:
476      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
477    }    }
478    return result.borrowReadyPtr();    return result.borrowReadyPtr();
479  }  }
# Line 480  DataLazy::collapse() Line 500  DataLazy::collapse()
500  }  }
501    
502  /*  /*
503    \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.
504    \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.
505    \param v A vector to store intermediate results.    \param v A vector to store intermediate results.
506    \param offset Index in v to begin storing results.    \param offset Index in v to begin storing results.
# Line 508  DataLazy::resolveUnary(ValueType& v, siz Line 528  DataLazy::resolveUnary(ValueType& v, siz
528    switch (m_op)    switch (m_op)
529    {    {
530      case SIN:        case SIN:  
531      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
532      break;      break;
533      case COS:      case COS:
534      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
535      break;      break;
536      case TAN:      case TAN:
537      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
538      break;      break;
539      case ASIN:      case ASIN:
540      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
541      break;      break;
542      case ACOS:      case ACOS:
543      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
544      break;      break;
545      case ATAN:      case ATAN:
546      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
547      break;      break;
548      case SINH:      case SINH:
549      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
550      break;      break;
551      case COSH:      case COSH:
552      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
553      break;      break;
554      case TANH:      case TANH:
555      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
556      break;      break;
557      case ERF:      case ERF:
558  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
559      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
560  #else  #else
561      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
562      break;      break;
563  #endif  #endif
564     case ASINH:     case ASINH:
565  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
566      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
567  #else  #else
568      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
569  #endif    #endif  
570      break;      break;
571     case ACOSH:     case ACOSH:
572  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
573      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
574  #else  #else
575      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
576  #endif    #endif  
577      break;      break;
578     case ATANH:     case ATANH:
579  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
580      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
581  #else  #else
582      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
583  #endif    #endif  
584      break;      break;
585      case LOG10:      case LOG10:
586      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
587      break;      break;
588      case LOG:      case LOG:
589      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
590      break;      break;
591      case SIGN:      case SIGN:
592      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
593      break;      break;
594      case ABS:      case ABS:
595      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
596      break;      break;
597      case NEG:      case NEG:
598      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 583  DataLazy::resolveUnary(ValueType& v, siz Line 603  DataLazy::resolveUnary(ValueType& v, siz
603      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
604      break;      break;
605      case EXP:      case EXP:
606      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
607      break;      break;
608      case SQRT:      case SQRT:
609      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
610      break;      break;
611      case RECIP:      case RECIP:
612      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 611  DataLazy::resolveUnary(ValueType& v, siz Line 631  DataLazy::resolveUnary(ValueType& v, siz
631  }  }
632    
633    
634    /*
635      \brief Compute the value of the expression (unary operation) for the given sample.
636      \return Vector which stores the value of the subexpression for the given sample.
637      \param v A vector to store intermediate results.
638      \param offset Index in v to begin storing results.
639      \param sampleNo Sample number to evaluate.
640      \param roffset (output parameter) the offset in the return vector where the result begins.
641    
642      The return value will be an existing vector so do not deallocate it.
643      If the result is stored in v it should be stored at the offset given.
644      Everything from offset to the end of v should be considered available for this method to use.
645    */
646    DataTypes::ValueType*
647    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
648    {
649        // we assume that any collapsing has been done before we get here
650        // since we only have one argument we don't need to think about only
651        // processing single points.
652      if (m_readytype!='E')
653      {
654        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
655      }
656        // since we can't write the result over the input, we need a result offset further along
657      size_t subroffset=roffset+m_samplesize;
658      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
659      roffset=offset;
660      switch (m_op)
661      {
662        case SYM:
663        DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
664        break;
665        case NSYM:
666        DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
667        break;
668        default:
669        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
670      }
671      return &v;
672    }
673    
674    
675    
676    
677  #define PROC_OP(X) \  #define PROC_OP(TYPE,X)                               \
678      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int i=0;i<steps;++i,resultp+=resultStep) \
679      { \      { \
680         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
681         lroffset+=leftStep; \         lroffset+=leftStep; \
682         rroffset+=rightStep; \         rroffset+=rightStep; \
683      }      }
# Line 672  cout << "Resolve binary: " << toString() Line 732  cout << "Resolve binary: " << toString()
732    switch(m_op)    switch(m_op)
733    {    {
734      case ADD:      case ADD:
735      PROC_OP(plus<double>());          PROC_OP(NO_ARG,plus<double>());
736      break;      break;
737      case SUB:      case SUB:
738      PROC_OP(minus<double>());      PROC_OP(NO_ARG,minus<double>());
739      break;      break;
740      case MUL:      case MUL:
741      PROC_OP(multiplies<double>());      PROC_OP(NO_ARG,multiplies<double>());
742      break;      break;
743      case DIV:      case DIV:
744      PROC_OP(divides<double>());      PROC_OP(NO_ARG,divides<double>());
745      break;      break;
746      case POW:      case POW:
747      PROC_OP(::pow);         PROC_OP(double (double,double),::pow);
748      break;      break;
749      default:      default:
750      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
# Line 843  DataLazy::intoString(ostringstream& oss) Line 903  DataLazy::intoString(ostringstream& oss)
903      oss << ')';      oss << ')';
904      break;      break;
905    case G_UNARY:    case G_UNARY:
906      case G_NP1OUT:
907      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
908      m_left->intoString(oss);      m_left->intoString(oss);
909      oss << ')';      oss << ')';

Legend:
Removed from v.1943  
changed lines
  Added in v.2050

  ViewVC Help
Powered by ViewVC 1.1.26