/[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 2199 by jfenwick, Thu Jan 8 06:10:52 2009 UTC revision 2496 by jfenwick, Fri Jun 26 06:09:47 2009 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    #include "EscriptParams.h"
32    
33  #include <iomanip>      // for some fancy formatting in debug  #include <iomanip>      // for some fancy formatting in debug
34    
35  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
# Line 40  bool privdebug=false; Line 42  bool privdebug=false;
42  #define DISABLEDEBUG privdebug=false;  #define DISABLEDEBUG privdebug=false;
43  }  }
44    
45  #define SIZELIMIT  // #define SIZELIMIT
46  // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}  // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
47    //#define SIZELIMIT if ((m_height>7) || (m_children>127)) {resolveToIdentity();}
48    
49    #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}
50    
51  /*  /*
52  How does DataLazy work?  How does DataLazy work?
# Line 115  enum ES_opgroup Line 119  enum ES_opgroup
119     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
120     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
121     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
122     G_TENSORPROD     // general tensor product     G_TENSORPROD,    // general tensor product
123       G_NP1OUT_2P      // non-pointwise op with one output requiring two params
124  };  };
125    
126    
# Line 129  string ES_opstrings[]={"UNKNOWN","IDENTI Line 134  string ES_opstrings[]={"UNKNOWN","IDENTI
134              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
135              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
136              "prod",              "prod",
137              "transpose", "trace"};              "transpose", "trace",
138  int ES_opcount=40;              "swapaxes"};
139    int ES_opcount=41;
140  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,
141              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
142              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 139  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT Line 145  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT
145              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35
146              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
147              G_TENSORPROD,              G_TENSORPROD,
148              G_NP1OUT_P, G_NP1OUT_P};              G_NP1OUT_P, G_NP1OUT_P,
149                G_NP1OUT_2P};
150  inline  inline
151  ES_opgroup  ES_opgroup
152  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 281  resultShape(DataAbstract_ptr left, ES_op Line 288  resultShape(DataAbstract_ptr left, ES_op
288      }      }
289  }  }
290    
291    DataTypes::ShapeType
292    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
293    {
294         // This code taken from the Data.cpp swapaxes() method
295         // Some of the checks are probably redundant here
296         int axis0_tmp,axis1_tmp;
297         const DataTypes::ShapeType& s=left->getShape();
298         DataTypes::ShapeType out_shape;
299         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
300         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
301         int rank=left->getRank();
302         if (rank<2) {
303            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
304         }
305         if (axis0<0 || axis0>rank-1) {
306            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
307         }
308         if (axis1<0 || axis1>rank-1) {
309             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
310         }
311         if (axis0 == axis1) {
312             throw DataException("Error - Data::swapaxes: axis indices must be different.");
313         }
314         if (axis0 > axis1) {
315             axis0_tmp=axis1;
316             axis1_tmp=axis0;
317         } else {
318             axis0_tmp=axis0;
319             axis1_tmp=axis1;
320         }
321         for (int i=0; i<rank; i++) {
322           if (i == axis0_tmp) {
323              out_shape.push_back(s[axis1_tmp]);
324           } else if (i == axis1_tmp) {
325              out_shape.push_back(s[axis0_tmp]);
326           } else {
327              out_shape.push_back(s[i]);
328           }
329         }
330        return out_shape;
331    }
332    
333    
334  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
335  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
336  // 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 348  GTPShape(DataAbstract_ptr left, DataAbst Line 398  GTPShape(DataAbstract_ptr left, DataAbst
398  // 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
399  // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.  // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
400  // The same goes for G_TENSORPROD  // The same goes for G_TENSORPROD
401  // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefths.  // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
402  // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined  // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
403  // multiple times  // multiple times
404  int  int
# Line 363  calcBuffs(const DataLazy_ptr& left, cons Line 413  calcBuffs(const DataLazy_ptr& left, cons
413     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
414     case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
415     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
416       case G_NP1OUT_2P: return 1+max(left->getBuffsRequired(),1);
417     default:     default:
418      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
419     }     }
# Line 397  DataLazy::DataLazy(DataAbstract_ptr p) Line 448  DataLazy::DataLazy(DataAbstract_ptr p)
448     }     }
449     else     else
450     {     {
451        p->makeLazyShared();
452      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
453      makeIdentity(dr);      makeIdentity(dr);
454  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
# Line 639  DataLazy::DataLazy(DataAbstract_ptr left Line 691  DataLazy::DataLazy(DataAbstract_ptr left
691  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
692  }  }
693    
694    
695    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
696        : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
697        m_op(op),
698        m_axis_offset(axis0),
699        m_transpose(axis1),
700        m_tol(0)
701    {
702       if ((getOpgroup(op)!=G_NP1OUT_2P))
703       {
704        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
705       }
706       DataLazy_ptr lleft;
707       if (!left->isLazy())
708       {
709        lleft=DataLazy_ptr(new DataLazy(left));
710       }
711       else
712       {
713        lleft=dynamic_pointer_cast<DataLazy>(left);
714       }
715       m_readytype=lleft->m_readytype;
716       m_left=lleft;
717       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
718       m_samplesize=getNumDPPSample()*getNoValues();
719       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
720       m_children=m_left->m_children+1;
721       m_height=m_left->m_height+1;
722       SIZELIMIT
723    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
724    }
725    
726  DataLazy::~DataLazy()  DataLazy::~DataLazy()
727  {  {
728  }  }
# Line 657  DataLazy::getMaxSampleSize() const Line 741  DataLazy::getMaxSampleSize() const
741      return m_maxsamplesize;      return m_maxsamplesize;
742  }  }
743    
744    
745    
746    size_t
747    DataLazy::getSampleBufferSize() const
748    {
749        return m_maxsamplesize*(max(1,m_buffsRequired));
750    }
751    
752  /*  /*
753    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
754    This does the work for the collapse method.    This does the work for the collapse method.
# Line 796  DataLazy::collapseToReady() Line 888  DataLazy::collapseToReady()
888      case TRACE:      case TRACE:
889      result=left.trace(m_axis_offset);      result=left.trace(m_axis_offset);
890      break;      break;
891        case SWAP:
892        result=left.swapaxes(m_axis_offset, m_transpose);
893        break;
894      default:      default:
895      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)+".");
896    }    }
# Line 1091  LAZYDEBUG(cerr << "Result of trace=" << Line 1186  LAZYDEBUG(cerr << "Result of trace=" <<
1186  }  }
1187    
1188    
1189    /*
1190      \brief Compute the value of the expression (unary operation with int params) for the given sample.
1191      \return Vector which stores the value of the subexpression for the given sample.
1192      \param v A vector to store intermediate results.
1193      \param offset Index in v to begin storing results.
1194      \param sampleNo Sample number to evaluate.
1195      \param roffset (output parameter) the offset in the return vector where the result begins.
1196    
1197      The return value will be an existing vector so do not deallocate it.
1198      If the result is stored in v it should be stored at the offset given.
1199      Everything from offset to the end of v should be considered available for this method to use.
1200    */
1201    DataTypes::ValueType*
1202    DataLazy::resolveNP1OUT_2P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1203    {
1204        // we assume that any collapsing has been done before we get here
1205        // since we only have one argument we don't need to think about only
1206        // processing single points.
1207      if (m_readytype!='E')
1208      {
1209        throw DataException("Programmer error - resolveNP1OUT_2P should only be called on expanded Data.");
1210      }
1211        // since we can't write the result over the input, we need a result offset further along
1212      size_t subroffset;
1213      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1214    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1215    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1216      roffset=offset;
1217      size_t loop=0;
1218      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1219      size_t outstep=getNoValues();
1220      size_t instep=m_left->getNoValues();
1221    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1222      switch (m_op)
1223      {
1224        case SWAP:
1225        for (loop=0;loop<numsteps;++loop)
1226        {
1227                DataMaths::swapaxes(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset, m_transpose);
1228            subroffset+=instep;
1229            offset+=outstep;
1230        }
1231        break;
1232        default:
1233        throw DataException("Programmer error - resolveNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1234      }
1235      return &v;
1236    }
1237    
1238    
1239    
1240  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1241      for (int j=0;j<onumsteps;++j)\      for (int j=0;j<onumsteps;++j)\
1242      {\      {\
# Line 1425  LAZYDEBUG(cout << "Resolve sample " << t Line 1571  LAZYDEBUG(cout << "Resolve sample " << t
1571    }    }
1572    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
1573    {    {
1574      const ValueType& vec=m_id->getVector();      const ValueType& vec=m_id->getVectorRO();
1575      if (m_readytype=='C')      if (m_readytype=='C')
1576      {      {
1577      roffset=0;      roffset=0;
# Line 1448  LAZYDEBUG(cout << "Finish  sample " << t Line 1594  LAZYDEBUG(cout << "Finish  sample " << t
1594    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1595    case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);    case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1596    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1597      case G_NP1OUT_2P: return resolveNP1OUT_2P(v, offset, sampleNo, roffset);
1598    default:    default:
1599      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)+".");
1600    }    }
1601    
1602  }  }
1603    
1604    const DataTypes::ValueType*
1605    DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
1606    {
1607    #ifdef _OPENMP
1608        int tid=omp_get_thread_num();
1609    #else
1610        int tid=0;
1611    #endif
1612        return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
1613    }
1614    
1615    
1616  // This needs to do the work of the idenity constructor  // This needs to do the work of the idenity constructor
1617  void  void
1618  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
# Line 1491  DataLazy::resolve() Line 1650  DataLazy::resolve()
1650    
1651  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1652  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
   
1653    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
1654    {    {
1655      collapse();      collapse();
# Line 1505  LAZYDEBUG(cout << "Buffers=" << m_buffsR Line 1663  LAZYDEBUG(cout << "Buffers=" << m_buffsR
1663      // storage to evaluate its expression      // storage to evaluate its expression
1664    int numthreads=1;    int numthreads=1;
1665  #ifdef _OPENMP  #ifdef _OPENMP
1666    numthreads=getNumberOfThreads();    numthreads=omp_get_max_threads();
1667  #endif  #endif
1668    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1669  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1670    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1671    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
1672    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1673    int sample;    int sample;
1674    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
# Line 1559  DataLazy::toString() const Line 1717  DataLazy::toString() const
1717  void  void
1718  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1719  {  {
1720  //   oss << "[" << m_children <<";"<<m_height <<"]";  //    oss << "[" << m_children <<";"<<m_height <<"]";
1721    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1722    {    {
1723    case G_IDENTITY:    case G_IDENTITY:
# Line 1603  DataLazy::intoString(ostringstream& oss) Line 1761  DataLazy::intoString(ostringstream& oss)
1761      m_right->intoString(oss);      m_right->intoString(oss);
1762      oss << ')';      oss << ')';
1763      break;      break;
1764      case G_NP1OUT_2P:
1765        oss << opToString(m_op) << '(';
1766        m_left->intoString(oss);
1767        oss << ", " << m_axis_offset << ", " << m_transpose;
1768        oss << ')';
1769        break;
1770    default:    default:
1771      oss << "UNKNOWN";      oss << "UNKNOWN";
1772    }    }
# Line 1699  DataLazy::getPointOffset(int sampleNo, Line 1863  DataLazy::getPointOffset(int sampleNo,
1863    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1864  }  }
1865    
1866  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
1867  // to zero, all the tags from all the DataTags would be in the result.  // I have decided to let Data:: handle this issue.
 // However since they all have the same value (0) whether they are there or not should not matter.  
 // So I have decided that for all types this method will create a constant 0.  
 // It can be promoted up as required.  
 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management  
 // but we can deal with that if it arrises.  
1868  void  void
1869  DataLazy::setToZero()  DataLazy::setToZero()
1870  {  {
1871    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
1872    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1873    m_op=IDENTITY;  //   m_op=IDENTITY;
1874    m_right.reset();    //   m_right.reset();  
1875    m_left.reset();  //   m_left.reset();
1876    m_readytype='C';  //   m_readytype='C';
1877    m_buffsRequired=1;  //   m_buffsRequired=1;
1878    
1879      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1880    
1881    }
1882    
1883    bool
1884    DataLazy::actsExpanded() const
1885    {
1886        return (m_readytype=='E');
1887  }  }
1888    
1889  }   // end namespace  }   // end namespace

Legend:
Removed from v.2199  
changed lines
  Added in v.2496

  ViewVC Help
Powered by ViewVC 1.1.26