/[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

revision 2082 by caltinay, Fri Nov 21 01:46:05 2008 UTC revision 2092 by jfenwick, Tue Nov 25 04:18:17 2008 UTC
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31    // #define LAZYDEBUG(X) X;
32    #define LAZYDEBUG(X)
33    
34  /*  /*
35  How does DataLazy work?  How does DataLazy work?
36  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 97  enum ES_opgroup Line 100  enum ES_opgroup
100     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
101     G_UNARY,     // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
102     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
103       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
104     G_TENSORPROD     // general tensor product     G_TENSORPROD     // general tensor product
105  };  };
106    
# Line 110  string ES_opstrings[]={"UNKNOWN","IDENTI Line 114  string ES_opstrings[]={"UNKNOWN","IDENTI
114              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
115              "1/","where>0","where<0","where>=0","where<=0",              "1/","where>0","where<0","where>=0","where<=0",
116              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
117              "prod"};              "prod",
118  int ES_opcount=36;              "transpose",
119                "trace"};
120    int ES_opcount=38;
121  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,
122              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
123              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
# Line 119  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT Line 125  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT
125              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
126              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,            // 33              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,            // 33
127              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
128              G_TENSORPROD};              G_TENSORPROD,
129                G_NP1OUT_P, G_NP1OUT_P};
130  inline  inline
131  ES_opgroup  ES_opgroup
132  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 177  resultShape(DataAbstract_ptr left, DataA Line 184  resultShape(DataAbstract_ptr left, DataA
184      return left->getShape();      return left->getShape();
185  }  }
186    
187    // return the shape for "op left"
188    
189    DataTypes::ShapeType
190    resultShape(DataAbstract_ptr left, ES_optype op)
191    {
192        switch(op)
193        {
194            case TRANS:
195            return left->getShape();
196        break;
197        case TRACE:
198            return DataTypes::scalarShape;
199        break;
200            default:
201        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
202        }
203    }
204    
205  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
206  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
207  // the majority of this code is copy pasted from C_General_Tensor_Product  // the majority of this code is copy pasted from C_General_Tensor_Product
# Line 197  GTPShape(DataAbstract_ptr left, DataAbst Line 222  GTPShape(DataAbstract_ptr left, DataAbst
222    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)  { start1 = rank1-axis_offset; }
223    else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }    else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
224    
225      if (rank0<axis_offset)
226      {
227        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
228      }
229    
230    // Adjust the shapes for transpose    // Adjust the shapes for transpose
231    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
# Line 226  GTPShape(DataAbstract_ptr left, DataAbst Line 255  GTPShape(DataAbstract_ptr left, DataAbst
255       for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z       for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
256       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
257    }    }
258    
259      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
260      {
261         ostringstream os;
262         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
263         throw DataException(os.str());
264      }
265    
266    return shape2;    return shape2;
267  }  }
268    
# Line 257  calcBuffs(const DataLazy_ptr& left, cons Line 294  calcBuffs(const DataLazy_ptr& left, cons
294     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
295     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY: return max(left->getBuffsRequired(),1);
296     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
297       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
298     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
299     default:     default:
300      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
# Line 305  DataLazy::DataLazy(DataAbstract_ptr p) Line 343  DataLazy::DataLazy(DataAbstract_ptr p)
343     m_buffsRequired=1;     m_buffsRequired=1;
344     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
345     m_maxsamplesize=m_samplesize;     m_maxsamplesize=m_samplesize;
346  cout << "(1)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
347  }  }
348    
349    
# Line 398  DataLazy::DataLazy(DataAbstract_ptr left Line 436  DataLazy::DataLazy(DataAbstract_ptr left
436     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
437     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
438     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
439  cout << "(3)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
440  }  }
441    
442  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
# Line 462  DataLazy::DataLazy(DataAbstract_ptr left Line 500  DataLazy::DataLazy(DataAbstract_ptr left
500     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
501     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
502     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
503  cout << "(4)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
504    }
505    
506    
507    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
508        : parent(left->getFunctionSpace(), resultShape(left,op)),
509        m_op(op),
510        m_axis_offset(axis_offset),
511        m_transpose(0)
512    {
513       if ((getOpgroup(op)!=G_NP1OUT_P))
514       {
515        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
516       }
517       DataLazy_ptr lleft;
518       if (!left->isLazy())
519       {
520        lleft=DataLazy_ptr(new DataLazy(left));
521       }
522       else
523       {
524        lleft=dynamic_pointer_cast<DataLazy>(left);
525       }
526       m_readytype=lleft->m_readytype;
527       m_left=lleft;
528       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
529       m_samplesize=getNumDPPSample()*getNoValues();
530       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
531    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
532  }  }
533    
534    
# Line 611  DataLazy::collapseToReady() Line 677  DataLazy::collapseToReady()
677      case PROD:      case PROD:
678      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
679      break;      break;
680        case TRANS:
681        result=left.transpose(m_axis_offset);
682        break;
683        case TRACE:
684        result=left.trace(m_axis_offset);
685        break;
686      default:      default:
687      throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
688    }    }
# Line 810  DataLazy::resolveNP1OUT(ValueType& v, si Line 882  DataLazy::resolveNP1OUT(ValueType& v, si
882    return &v;    return &v;
883  }  }
884    
885    /*
886      \brief Compute the value of the expression (unary operation) for the given sample.
887      \return Vector which stores the value of the subexpression for the given sample.
888      \param v A vector to store intermediate results.
889      \param offset Index in v to begin storing results.
890      \param sampleNo Sample number to evaluate.
891      \param roffset (output parameter) the offset in the return vector where the result begins.
892    
893      The return value will be an existing vector so do not deallocate it.
894      If the result is stored in v it should be stored at the offset given.
895      Everything from offset to the end of v should be considered available for this method to use.
896    */
897    DataTypes::ValueType*
898    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
899    {
900        // we assume that any collapsing has been done before we get here
901        // since we only have one argument we don't need to think about only
902        // processing single points.
903      if (m_readytype!='E')
904      {
905        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
906      }
907        // since we can't write the result over the input, we need a result offset further along
908      size_t subroffset=roffset+m_samplesize;
909      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
910      roffset=offset;
911      switch (m_op)
912      {
913        case TRACE:
914             DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
915        break;
916        case TRANS:
917             DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
918        break;
919        default:
920        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
921      }
922      return &v;
923    }
924    
925    
926  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
# Line 845  DataLazy::resolveNP1OUT(ValueType& v, si Line 955  DataLazy::resolveNP1OUT(ValueType& v, si
955  DataTypes::ValueType*  DataTypes::ValueType*
956  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
957  {  {
958  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
959    
960    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
961      // first work out which of the children are expanded      // first work out which of the children are expanded
962    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
963    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
964      if (!leftExp && !rightExp)
965      {
966        throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
967      }
968      bool leftScalar=(m_left->getRank()==0);
969      bool rightScalar=(m_right->getRank()==0);
970    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
971    int steps=(bigloops?1:getNumDPPSample());    int steps=(bigloops?1:getNumDPPSample());
972    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
973    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
974    {    {
975      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");      if (!leftScalar && !rightScalar)
976        {
977           throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
978        }
979      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
980      chunksize=1;    // for scalar      chunksize=1;    // for scalar
981    }        }    
982    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    int leftStep=((leftExp && (!rightExp || rightScalar))? m_right->getNoValues() : 0);
983    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=((rightExp && (!leftExp || leftScalar))? m_left->getNoValues() : 0);
984    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
985      // Get the values of sub-expressions      // Get the values of sub-expressions
986    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
# Line 911  cout << "Resolve binary: " << toString() Line 1030  cout << "Resolve binary: " << toString()
1030  DataTypes::ValueType*  DataTypes::ValueType*
1031  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1032  {  {
1033  cout << "Resolve TensorProduct: " << toString() << endl;  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1034    
1035    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
1036      // first work out which of the children are expanded      // first work out which of the children are expanded
# Line 965  cout << "Resolve TensorProduct: " << toS Line 1084  cout << "Resolve TensorProduct: " << toS
1084  const DataTypes::ValueType*  const DataTypes::ValueType*
1085  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1086  {  {
1087  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1088      // 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
1089    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1090    {    {
# Line 991  cout << "Resolve sample " << toString() Line 1110  cout << "Resolve sample " << toString()
1110    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
1111    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1112    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1113      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1114    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1115    default:    default:
1116      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)+".");
# Line 1004  DataReady_ptr Line 1124  DataReady_ptr
1124  DataLazy::resolve()  DataLazy::resolve()
1125  {  {
1126    
1127  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1128  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1129    
1130    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
1131    {    {
# Line 1023  cout << "Buffers=" << m_buffsRequired << Line 1143  cout << "Buffers=" << m_buffsRequired <<
1143    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
1144  #endif  #endif
1145    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1146  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1147    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1148    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1149    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
# Line 1035  cout << "Buffer created with size=" << v Line 1155  cout << "Buffer created with size=" << v
1155    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1156    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1157    {    {
1158  cout << "################################# " << sample << endl;  LAZYDEBUG(cout << "################################# " << sample << endl;)
1159  #ifdef _OPENMP  #ifdef _OPENMP
1160      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1161  #else  #else
1162      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.
1163  #endif  #endif
1164  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1165      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1166  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)
1167      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
1168      {      {
1169      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1170      }      }
1171  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << "*********************************" << endl;)
1172    }    }
1173    return resptr;    return resptr;
1174  }  }
# Line 1096  DataLazy::intoString(ostringstream& oss) Line 1216  DataLazy::intoString(ostringstream& oss)
1216      break;      break;
1217    case G_UNARY:    case G_UNARY:
1218    case G_NP1OUT:    case G_NP1OUT:
1219      case G_NP1OUT_P:
1220      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1221      m_left->intoString(oss);      m_left->intoString(oss);
1222      oss << ')';      oss << ')';

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

  ViewVC Help
Powered by ViewVC 1.1.26