/[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 2147 by jfenwick, Wed Dec 10 04:41:26 2008 UTC revision 2177 by jfenwick, Wed Dec 17 23:51:23 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;  // #define LAZYDEBUG(X) if (privdebug){X;}
32  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
33    namespace
34    {
35    bool privdebug=false;
36    
37    #define ENABLEDEBUG privdebug=true;
38    #define DISABLEDEBUG privdebug=false;
39    }
40    
41    #define SIZELIMIT
42    // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
43    
44    
45  /*  /*
46  How does DataLazy work?  How does DataLazy work?
# Line 266  GTPShape(DataAbstract_ptr left, DataAbst Line 277  GTPShape(DataAbstract_ptr left, DataAbst
277    return shape2;    return shape2;
278  }  }
279    
   
 // determine the number of points in the result of "left op right"  
 // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here  
 // size_t  
 // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
 // {  
 //    switch (getOpgroup(op))  
 //    {  
 //    case G_BINARY: return left->getLength();  
 //    case G_UNARY: return left->getLength();  
 //    case G_NP1OUT: return left->getLength();  
 //    default:  
 //  throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");  
 //    }  
 // }  
   
280  // 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
281  // 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.
282  // The same goes for G_TENSORPROD  // The same goes for G_TENSORPROD
283    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefths.
284    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
285    // multiple times
286  int  int
287  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
288  {  {
289     switch(getOpgroup(op))     switch(getOpgroup(op))
290     {     {
291     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
292     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
293     case G_UNARY:     case G_UNARY:
294     case G_UNARY_P: return max(left->getBuffsRequired(),1);     case G_UNARY_P: return max(left->getBuffsRequired(),1);
295     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
# Line 320  opToString(ES_optype op) Line 318  opToString(ES_optype op)
318    
319    
320  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
321      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
     m_op(IDENTITY),  
     m_axis_offset(0),  
     m_transpose(0),  
     m_SL(0), m_SM(0), m_SR(0)  
322  {  {
323     if (p->isLazy())     if (p->isLazy())
324     {     {
# Line 335  DataLazy::DataLazy(DataAbstract_ptr p) Line 329  DataLazy::DataLazy(DataAbstract_ptr p)
329     }     }
330     else     else
331     {     {
332      m_id=dynamic_pointer_cast<DataReady>(p);      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
333      if(p->isConstant()) {m_readytype='C';}      makeIdentity(dr);
     else if(p->isExpanded()) {m_readytype='E';}  
     else if (p->isTagged()) {m_readytype='T';}  
     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}  
334     }     }
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_maxsamplesize=m_samplesize;  
335  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
336  }  }
337    
# Line 376  DataLazy::DataLazy(DataAbstract_ptr left Line 364  DataLazy::DataLazy(DataAbstract_ptr left
364     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
365     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
366     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
367       m_children=m_left->m_children+1;
368       m_height=m_left->m_height+1;
369       SIZELIMIT
370  }  }
371    
372    
# Line 437  DataLazy::DataLazy(DataAbstract_ptr left Line 428  DataLazy::DataLazy(DataAbstract_ptr left
428     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
429     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
430     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
431       m_children=m_left->m_children+m_right->m_children+2;
432       m_height=max(m_left->m_height,m_right->m_height)+1;
433       SIZELIMIT
434  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
435  }  }
436    
# Line 501  DataLazy::DataLazy(DataAbstract_ptr left Line 495  DataLazy::DataLazy(DataAbstract_ptr left
495     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
496     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
497     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
498       m_children=m_left->m_children+m_right->m_children+2;
499       m_height=max(m_left->m_height,m_right->m_height)+1;
500       SIZELIMIT
501  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
502  }  }
503    
# Line 530  DataLazy::DataLazy(DataAbstract_ptr left Line 527  DataLazy::DataLazy(DataAbstract_ptr left
527     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
528     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
529     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
530       m_children=m_left->m_children+1;
531       m_height=m_left->m_height+1;
532       SIZELIMIT
533  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
534  }  }
535    
# Line 558  DataLazy::DataLazy(DataAbstract_ptr left Line 558  DataLazy::DataLazy(DataAbstract_ptr left
558     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
559     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
560     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
561       m_children=m_left->m_children+1;
562       m_height=m_left->m_height+1;
563       SIZELIMIT
564  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
565  }  }
566    
# Line 912  DataLazy::resolveNP1OUT(ValueType& v, si Line 915  DataLazy::resolveNP1OUT(ValueType& v, si
915    }    }
916      // since we can't write the result over the input, we need a result offset further along      // since we can't write the result over the input, we need a result offset further along
917    size_t subroffset=roffset+m_samplesize;    size_t subroffset=roffset+m_samplesize;
918    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);  LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
919      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
920    roffset=offset;    roffset=offset;
921      size_t loop=0;
922      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
923      size_t step=getNoValues();
924    switch (m_op)    switch (m_op)
925    {    {
926      case SYM:      case SYM:
927      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
928        {
929            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
930            subroffset+=step;
931            offset+=step;
932        }
933      break;      break;
934      case NSYM:      case NSYM:
935      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
936        {
937            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
938            subroffset+=step;
939            offset+=step;
940        }
941      break;      break;
942      default:      default:
943      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
# Line 951  DataLazy::resolveNP1OUT_P(ValueType& v, Line 968  DataLazy::resolveNP1OUT_P(ValueType& v,
968      throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");      throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
969    }    }
970      // since we can't write the result over the input, we need a result offset further along      // since we can't write the result over the input, we need a result offset further along
971    size_t subroffset=roffset+m_samplesize;    size_t subroffset;
972    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);    const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
973    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
974    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
975    roffset=offset;    roffset=offset;
976      size_t loop=0;
977      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
978      size_t outstep=getNoValues();
979      size_t instep=m_left->getNoValues();
980    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
981    switch (m_op)    switch (m_op)
982    {    {
983      case TRACE:      case TRACE:
984           DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
985        {
986    size_t zz=sampleNo*getNumDPPSample()+loop;
987    if (zz==5800)
988    {
989    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
990    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
991    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
992    LAZYDEBUG(cerr << subroffset << endl;)
993    LAZYDEBUG(cerr << "output=" << offset << endl;)
994    }
995                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
996    if (zz==5800)
997    {
998    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
999    }
1000            subroffset+=instep;
1001            offset+=outstep;
1002        }
1003      break;      break;
1004      case TRANS:      case TRANS:
1005           DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
1006        {
1007                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1008            subroffset+=instep;
1009            offset+=outstep;
1010        }
1011      break;      break;
1012      default:      default:
1013      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
# Line 970  DataLazy::resolveNP1OUT_P(ValueType& v, Line 1017  DataLazy::resolveNP1OUT_P(ValueType& v,
1017    
1018    
1019  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1020      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int j=0;j<onumsteps;++j)\
1021      { \      {\
1022         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1023         lroffset+=leftStep; \        { \
1024         rroffset+=rightStep; \  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1025    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1026             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1027    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1028             lroffset+=leftstep; \
1029             rroffset+=rightstep; \
1030          }\
1031          lroffset+=oleftstep;\
1032          rroffset+=orightstep;\
1033      }      }
1034    
1035  /*  /*
# Line 1013  LAZYDEBUG(cout << "Resolve binary: " << Line 1068  LAZYDEBUG(cout << "Resolve binary: " <<
1068    }    }
1069    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1070    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1071    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
   int steps=(bigloops?1:getNumDPPSample());  
   size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint  
   if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops  
1072    {    {
1073      if (!leftScalar && !rightScalar)      throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1074      }
1075      size_t leftsize=m_left->getNoValues();
1076      size_t rightsize=m_right->getNoValues();
1077      size_t chunksize=1;           // how many doubles will be processed in one go
1078      int leftstep=0;       // how far should the left offset advance after each step
1079      int rightstep=0;
1080      int numsteps=0;       // total number of steps for the inner loop
1081      int oleftstep=0;  // the o variables refer to the outer loop
1082      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1083      int onumsteps=1;
1084      
1085      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1086      bool RES=(rightExp && rightScalar);
1087      bool LS=(!leftExp && leftScalar); // left is a single scalar
1088      bool RS=(!rightExp && rightScalar);
1089      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1090      bool RN=(!rightExp && !rightScalar);
1091      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1092      bool REN=(rightExp && !rightScalar);
1093    
1094      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1095      {
1096        chunksize=m_left->getNumDPPSample()*leftsize;
1097        leftstep=0;
1098        rightstep=0;
1099        numsteps=1;
1100      }
1101      else if (LES || RES)
1102      {
1103        chunksize=1;
1104        if (LES)        // left is an expanded scalar
1105      {      {
1106         throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");          if (RS)
1107            {
1108               leftstep=1;
1109               rightstep=0;
1110               numsteps=m_left->getNumDPPSample();
1111            }
1112            else        // RN or REN
1113            {
1114               leftstep=0;
1115               oleftstep=1;
1116               rightstep=1;
1117               orightstep=(RN ? -(int)rightsize : 0);
1118               numsteps=rightsize;
1119               onumsteps=m_left->getNumDPPSample();
1120            }
1121      }      }
1122      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());      else        // right is an expanded scalar
1123      chunksize=1;    // for scalar      {
1124    }              if (LS)
1125    int leftStep=((leftExp && (!rightExp || rightScalar))? m_right->getNoValues() : 0);          {
1126    int rightStep=((rightExp && (!leftExp || leftScalar))? m_left->getNoValues() : 0);             rightstep=1;
1127    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0             leftstep=0;
1128               numsteps=m_right->getNumDPPSample();
1129            }
1130            else
1131            {
1132               rightstep=0;
1133               orightstep=1;
1134               leftstep=1;
1135               oleftstep=(LN ? -(int)leftsize : 0);
1136               numsteps=leftsize;
1137               onumsteps=m_right->getNumDPPSample();
1138            }
1139        }
1140      }
1141      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1142      {
1143        if (LEN)    // and Right will be a single value
1144        {
1145            chunksize=rightsize;
1146            leftstep=rightsize;
1147            rightstep=0;
1148            numsteps=m_left->getNumDPPSample();
1149            if (RS)
1150            {
1151               numsteps*=leftsize;
1152            }
1153        }
1154        else    // REN
1155        {
1156            chunksize=leftsize;
1157            rightstep=leftsize;
1158            leftstep=0;
1159            numsteps=m_right->getNumDPPSample();
1160            if (LS)
1161            {
1162               numsteps*=rightsize;
1163            }
1164        }
1165      }
1166    
1167      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1168      // Get the values of sub-expressions      // Get the values of sub-expressions
1169    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1170    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note      // calcBufss for why we can't put offset as the 2nd param above
1171      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1172      // the right child starts further along.      // the right child starts further along.
1173    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1174    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1175    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1176    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1177    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1178    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1179    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1180    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
1181    switch(m_op)    switch(m_op)
1182    {    {
# Line 1058  LAZYDEBUG(cout << "Resolve binary: " << Line 1203  LAZYDEBUG(cout << "Resolve binary: " <<
1203  }  }
1204    
1205    
1206    
1207  /*  /*
1208    \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.
1209    \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 1087  LAZYDEBUG(cout << "Resolve TensorProduct Line 1233  LAZYDEBUG(cout << "Resolve TensorProduct
1233    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
1234    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
1235      // Get the values of sub-expressions (leave a gap of one sample for the result).      // Get the values of sub-expressions (leave a gap of one sample for the result).
1236    const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);    int gap=offset+m_left->getMaxSampleSize();    // actually only needs to be m_left->m_samplesize
1237    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);    const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1238      gap+=m_right->getMaxSampleSize();
1239      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1240    LAZYDEBUG(cout << "Post sub calls: " << toString() << endl;)
1241    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1242    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1243    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1244    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1245    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1246    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
1247    switch(m_op)    switch(m_op)
1248    {    {
1249      case PROD:      case PROD:
1250      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1251      {      {
1252    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1253    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1254    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1255            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1256            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1257    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1258    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1259            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1260        lroffset+=leftStep;        lroffset+=leftStep;
1261        rroffset+=rightStep;        rroffset+=rightStep;
# Line 1142  LAZYDEBUG(cout << "Resolve sample " << t Line 1301  LAZYDEBUG(cout << "Resolve sample " << t
1301      if (m_readytype=='C')      if (m_readytype=='C')
1302      {      {
1303      roffset=0;      roffset=0;
1304    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1305      return &(vec);      return &(vec);
1306      }      }
1307      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1308    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1309      return &(vec);      return &(vec);
1310    }    }
1311    if (m_readytype!='E')    if (m_readytype!='E')
# Line 1162  LAZYDEBUG(cout << "Resolve sample " << t Line 1323  LAZYDEBUG(cout << "Resolve sample " << t
1323    default:    default:
1324      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)+".");
1325    }    }
1326    
1327  }  }
1328    
1329    // This needs to do the work of the idenity constructor
1330    void
1331    DataLazy::resolveToIdentity()
1332    {
1333       if (m_op==IDENTITY)
1334        return;
1335       DataReady_ptr p=resolve();
1336       makeIdentity(p);
1337    }
1338    
1339    void DataLazy::makeIdentity(const DataReady_ptr& p)
1340    {
1341       m_op=IDENTITY;
1342       m_axis_offset=0;
1343       m_transpose=0;
1344       m_SL=m_SM=m_SR=0;
1345       m_children=m_height=0;
1346       m_id=p;
1347       if(p->isConstant()) {m_readytype='C';}
1348       else if(p->isExpanded()) {m_readytype='E';}
1349       else if (p->isTagged()) {m_readytype='T';}
1350       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1351       m_buffsRequired=1;
1352       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1353       m_maxsamplesize=m_samplesize;
1354       m_left.reset();
1355       m_right.reset();
1356    }
1357    
1358  // To simplify the memory management, all threads operate on one large vector, rather than one each.  // To simplify the memory management, all threads operate on one large vector, rather than one each.
1359  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
# Line 1199  LAZYDEBUG(cout << "Buffer created with s Line 1389  LAZYDEBUG(cout << "Buffer created with s
1389    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1390    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1391    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1392    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1393    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1394    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1395    {    {
1396          if (sample==0)  {ENABLEDEBUG}
1397    
1398    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1399  LAZYDEBUG(cout << "################################# " << sample << endl;)  LAZYDEBUG(cout << "################################# " << sample << endl;)
1400  #ifdef _OPENMP  #ifdef _OPENMP
1401      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
# Line 1209  LAZYDEBUG(cout << "##################### Line 1403  LAZYDEBUG(cout << "#####################
1403      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.
1404  #endif  #endif
1405  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1406    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1407      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1408  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1409      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
1410      {      {
1411    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1412      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1413      }      }
1414    LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1415  LAZYDEBUG(cerr << "*********************************" << endl;)  LAZYDEBUG(cerr << "*********************************" << endl;)
1416        DISABLEDEBUG
1417    }    }
1418    return resptr;    return resptr;
1419  }  }
# Line 1233  DataLazy::toString() const Line 1431  DataLazy::toString() const
1431  void  void
1432  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1433  {  {
1434      oss << "[" << m_children <<";"<<m_height <<"]";
1435    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1436    {    {
1437    case G_IDENTITY:    case G_IDENTITY:

Legend:
Removed from v.2147  
changed lines
  Added in v.2177

  ViewVC Help
Powered by ViewVC 1.1.26