/[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 2152 by jfenwick, Thu Dec 11 04:26:42 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 70  The convention that I use, is that the r Line 73  The convention that I use, is that the r
73  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.
74  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.
75    
76  To add a new operator you need to do the following (plus anything I might have forgotten):  To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
77  1) Add to the ES_optype.  1) Add to the ES_optype.
78  2) determine what opgroup your operation belongs to (X)  2) determine what opgroup your operation belongs to (X)
79  3) add a string for the op to the end of ES_opstrings  3) add a string for the op to the end of ES_opstrings
# Line 96  enum ES_opgroup Line 99  enum ES_opgroup
99     G_IDENTITY,     G_IDENTITY,
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_UNARY_P,       // pointwise operations with one argument, requiring a parameter
103     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
104       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
105     G_TENSORPROD     // general tensor product     G_TENSORPROD     // general tensor product
106  };  };
107    
# Line 108  string ES_opstrings[]={"UNKNOWN","IDENTI Line 113  string ES_opstrings[]={"UNKNOWN","IDENTI
113              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
114              "asinh","acosh","atanh",              "asinh","acosh","atanh",
115              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
116              "1/","where>0","where<0","where>=0","where<=0",              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
117              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
118              "prod"};              "prod",
119  int ES_opcount=36;              "transpose", "trace"};
120    int ES_opcount=40;
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
124              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
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, G_UNARY_P, G_UNARY_P,      // 35
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 255  calcBuffs(const DataLazy_ptr& left, cons Line 292  calcBuffs(const DataLazy_ptr& left, cons
292     {     {
293     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
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:
296       case G_UNARY_P: return max(left->getBuffsRequired(),1);
297     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
298       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
299     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
300     default:     default:
301      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 344  DataLazy::DataLazy(DataAbstract_ptr p)
344     m_buffsRequired=1;     m_buffsRequired=1;
345     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
346     m_maxsamplesize=m_samplesize;     m_maxsamplesize=m_samplesize;
347  cout << "(1)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
348  }  }
349    
350    
# Line 398  DataLazy::DataLazy(DataAbstract_ptr left Line 437  DataLazy::DataLazy(DataAbstract_ptr left
437     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
438     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
439     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
440  cout << "(3)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
441  }  }
442    
443  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 501  DataLazy::DataLazy(DataAbstract_ptr left
501     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
502     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
503     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
504  cout << "(4)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
505    }
506    
507    
508    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
509        : parent(left->getFunctionSpace(), resultShape(left,op)),
510        m_op(op),
511        m_axis_offset(axis_offset),
512        m_transpose(0),
513        m_tol(0)
514    {
515       if ((getOpgroup(op)!=G_NP1OUT_P))
516       {
517        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
518       }
519       DataLazy_ptr lleft;
520       if (!left->isLazy())
521       {
522        lleft=DataLazy_ptr(new DataLazy(left));
523       }
524       else
525       {
526        lleft=dynamic_pointer_cast<DataLazy>(left);
527       }
528       m_readytype=lleft->m_readytype;
529       m_left=lleft;
530       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
531       m_samplesize=getNumDPPSample()*getNoValues();
532       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
533    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
534  }  }
535    
536    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
537        : parent(left->getFunctionSpace(), left->getShape()),
538        m_op(op),
539        m_axis_offset(0),
540        m_transpose(0),
541        m_tol(tol)
542    {
543       if ((getOpgroup(op)!=G_UNARY_P))
544       {
545        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
546       }
547       DataLazy_ptr lleft;
548       if (!left->isLazy())
549       {
550        lleft=DataLazy_ptr(new DataLazy(left));
551       }
552       else
553       {
554        lleft=dynamic_pointer_cast<DataLazy>(left);
555       }
556       m_readytype=lleft->m_readytype;
557       m_left=lleft;
558       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
559       m_samplesize=getNumDPPSample()*getNoValues();
560       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
561    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
562    }
563    
564  DataLazy::~DataLazy()  DataLazy::~DataLazy()
565  {  {
# Line 602  DataLazy::collapseToReady() Line 697  DataLazy::collapseToReady()
697      case LEZ:      case LEZ:
698      result=left.whereNonPositive();      result=left.whereNonPositive();
699      break;      break;
700        case NEZ:
701        result=left.whereNonZero(m_tol);
702        break;
703        case EZ:
704        result=left.whereZero(m_tol);
705        break;
706      case SYM:      case SYM:
707      result=left.symmetric();      result=left.symmetric();
708      break;      break;
# Line 611  DataLazy::collapseToReady() Line 712  DataLazy::collapseToReady()
712      case PROD:      case PROD:
713      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
714      break;      break;
715        case TRANS:
716        result=left.transpose(m_axis_offset);
717        break;
718        case TRACE:
719        result=left.trace(m_axis_offset);
720        break;
721      default:      default:
722      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)+".");
723    }    }
# Line 762  DataLazy::resolveUnary(ValueType& v, siz Line 869  DataLazy::resolveUnary(ValueType& v, siz
869      case LEZ:      case LEZ:
870      tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));      tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
871      break;      break;
872    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
873        case NEZ:
874        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
875        break;
876        case EZ:
877        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
878        break;
879    
880      default:      default:
881      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
# Line 770  DataLazy::resolveUnary(ValueType& v, siz Line 884  DataLazy::resolveUnary(ValueType& v, siz
884  }  }
885    
886    
887    
888    
889    
890    
891  /*  /*
892    \brief Compute the value of the expression (unary operation) for the given sample.    \brief Compute the value of the expression (unary operation) for the given sample.
893    \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.
# Line 810  DataLazy::resolveNP1OUT(ValueType& v, si Line 928  DataLazy::resolveNP1OUT(ValueType& v, si
928    return &v;    return &v;
929  }  }
930    
931    /*
932      \brief Compute the value of the expression (unary operation) for the given sample.
933      \return Vector which stores the value of the subexpression for the given sample.
934      \param v A vector to store intermediate results.
935      \param offset Index in v to begin storing results.
936      \param sampleNo Sample number to evaluate.
937      \param roffset (output parameter) the offset in the return vector where the result begins.
938    
939      The return value will be an existing vector so do not deallocate it.
940      If the result is stored in v it should be stored at the offset given.
941      Everything from offset to the end of v should be considered available for this method to use.
942    */
943    DataTypes::ValueType*
944    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
945    {
946        // we assume that any collapsing has been done before we get here
947        // since we only have one argument we don't need to think about only
948        // processing single points.
949      if (m_readytype!='E')
950      {
951        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
952      }
953        // since we can't write the result over the input, we need a result offset further along
954      size_t subroffset=roffset+m_samplesize;
955      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
956      roffset=offset;
957      switch (m_op)
958      {
959        case TRACE:
960             DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
961        break;
962        case TRANS:
963             DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
964        break;
965        default:
966        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
967      }
968      return &v;
969    }
970    
971    
972  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
973      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int j=0;j<onumsteps;++j)\
974      { \      {\
975         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
976         lroffset+=leftStep; \        { \
977         rroffset+=rightStep; \  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
978             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
979             lroffset+=leftstep; \
980             rroffset+=rightstep; \
981          }\
982          lroffset+=oleftstep;\
983          rroffset+=orightstep;\
984      }      }
985    
986  /*  /*
# Line 845  DataLazy::resolveNP1OUT(ValueType& v, si Line 1007  DataLazy::resolveNP1OUT(ValueType& v, si
1007  DataTypes::ValueType*  DataTypes::ValueType*
1008  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
1009  {  {
1010  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1011    
1012    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
1013      // first work out which of the children are expanded      // first work out which of the children are expanded
1014    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1015    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1016    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1017    int steps=(bigloops?1:getNumDPPSample());    {
1018    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint      throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1019    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1020    {    bool leftScalar=(m_left->getRank()==0);
1021      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1022      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1023      chunksize=1;    // for scalar    {
1024    }          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1025    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1026    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1027    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1028      size_t chunksize=1;           // how many doubles will be processed in one go
1029      int leftstep=0;       // how far should the left offset advance after each step
1030      int rightstep=0;
1031      int numsteps=0;       // total number of steps for the inner loop
1032      int oleftstep=0;  // the o variables refer to the outer loop
1033      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1034      int onumsteps=1;
1035      
1036      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1037      bool RES=(rightExp && rightScalar);
1038      bool LS=(!leftExp && leftScalar); // left is a single scalar
1039      bool RS=(!rightExp && rightScalar);
1040      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1041      bool RN=(!rightExp && !rightScalar);
1042      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1043      bool REN=(rightExp && !rightScalar);
1044    
1045      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1046      {
1047        chunksize=m_left->getNumDPPSample()*leftsize;
1048        leftstep=0;
1049        rightstep=0;
1050        numsteps=1;
1051      }
1052      else if (LES || RES)
1053      {
1054        chunksize=1;
1055        if (LES)        // left is an expanded scalar
1056        {
1057            if (RS)
1058            {
1059               leftstep=1;
1060               rightstep=0;
1061               numsteps=m_left->getNumDPPSample();
1062            }
1063            else        // RN or REN
1064            {
1065               leftstep=0;
1066               oleftstep=1;
1067               rightstep=1;
1068               orightstep=(RN?-rightsize:0);
1069               numsteps=rightsize;
1070               onumsteps=m_left->getNumDPPSample();
1071            }
1072        }
1073        else        // right is an expanded scalar
1074        {
1075            if (LS)
1076            {
1077               rightstep=1;
1078               leftstep=0;
1079               numsteps=m_right->getNumDPPSample();
1080            }
1081            else
1082            {
1083               rightstep=0;
1084               orightstep=1;
1085               leftstep=1;
1086               oleftstep=(LN?-leftsize:0);
1087               numsteps=leftsize;
1088               onumsteps=m_right->getNumDPPSample();
1089            }
1090        }
1091      }
1092      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1093      {
1094        if (LEN)    // and Right will be a single value
1095        {
1096            chunksize=rightsize;
1097            leftstep=rightsize;
1098            rightstep=0;
1099            numsteps=m_left->getNumDPPSample();
1100            if (RS)
1101            {
1102               numsteps*=leftsize;
1103            }
1104        }
1105        else    // REN
1106        {
1107            chunksize=leftsize;
1108            rightstep=leftsize;
1109            leftstep=0;
1110            numsteps=m_right->getNumDPPSample();
1111            if (LS)
1112            {
1113               numsteps*=rightsize;
1114            }
1115        }
1116      }
1117    
1118      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1119      // Get the values of sub-expressions      // Get the values of sub-expressions
1120    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
1121    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
1122      // the right child starts further along.      // the right child starts further along.
1123    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1124    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1125    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1126    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1127    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1128    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1129    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1130    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1131    switch(m_op)    switch(m_op)
1132    {    {
# Line 893  cout << "Resolve binary: " << toString() Line 1153  cout << "Resolve binary: " << toString()
1153  }  }
1154    
1155    
1156    
1157  /*  /*
1158    \brief Compute the value of the expression (tensor product) for the given sample.    \brief Compute the value of the expression (tensor product) for the given sample.
1159    \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.
# Line 911  cout << "Resolve binary: " << toString() Line 1172  cout << "Resolve binary: " << toString()
1172  DataTypes::ValueType*  DataTypes::ValueType*
1173  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
1174  {  {
1175  cout << "Resolve TensorProduct: " << toString() << endl;  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1176    
1177    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
1178      // 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 1226  cout << "Resolve TensorProduct: " << toS
1226  const DataTypes::ValueType*  const DataTypes::ValueType*
1227  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1228  {  {
1229  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1230      // 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
1231    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1232    {    {
# Line 977  cout << "Resolve sample " << toString() Line 1238  cout << "Resolve sample " << toString()
1238      if (m_readytype=='C')      if (m_readytype=='C')
1239      {      {
1240      roffset=0;      roffset=0;
1241    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1242      return &(vec);      return &(vec);
1243      }      }
1244      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1245    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1246      return &(vec);      return &(vec);
1247    }    }
1248    if (m_readytype!='E')    if (m_readytype!='E')
# Line 988  cout << "Resolve sample " << toString() Line 1251  cout << "Resolve sample " << toString()
1251    }    }
1252    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1253    {    {
1254    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY:
1255      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1256    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1257    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1258      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1259    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1260    default:    default:
1261      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)+".");
1262    }    }
1263    
1264  }  }
1265    
1266    
# Line 1004  DataReady_ptr Line 1270  DataReady_ptr
1270  DataLazy::resolve()  DataLazy::resolve()
1271  {  {
1272    
1273  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1274  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1275    
1276    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
1277    {    {
# Line 1023  cout << "Buffers=" << m_buffsRequired << Line 1289  cout << "Buffers=" << m_buffsRequired <<
1289    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
1290  #endif  #endif
1291    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1292  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1293    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1294    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1295    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
# Line 1032  cout << "Buffer created with size=" << v Line 1298  cout << "Buffer created with size=" << v
1298    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1299    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1300    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1301    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1302    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1303    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1304    {    {
1305  cout << "################################# " << sample << endl;  LAZYDEBUG(cout << "################################# " << sample << endl;)
1306  #ifdef _OPENMP  #ifdef _OPENMP
1307      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1308  #else  #else
1309      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.
1310  #endif  #endif
1311  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1312      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1313  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)
1314      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
1315      {      {
1316      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1317      }      }
1318  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << "*********************************" << endl;)
1319    }    }
1320    return resptr;    return resptr;
1321  }  }
# Line 1095  DataLazy::intoString(ostringstream& oss) Line 1362  DataLazy::intoString(ostringstream& oss)
1362      oss << ')';      oss << ')';
1363      break;      break;
1364    case G_UNARY:    case G_UNARY:
1365      case G_UNARY_P:
1366    case G_NP1OUT:    case G_NP1OUT:
1367      case G_NP1OUT_P:
1368      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1369      m_left->intoString(oss);      m_left->intoString(oss);
1370      oss << ')';      oss << ')';

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

  ViewVC Help
Powered by ViewVC 1.1.26