/[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 2195 by jfenwick, Wed Jan 7 04:13:52 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  // #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);
334      else if(p->isExpanded()) {m_readytype='E';}  cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;
     else if (p->isTagged()) {m_readytype='T';}  
     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}  
335     }     }
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_maxsamplesize=m_samplesize;  
336  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
337  }  }
338    
# Line 376  DataLazy::DataLazy(DataAbstract_ptr left Line 365  DataLazy::DataLazy(DataAbstract_ptr left
365     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
366     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
367     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
368       m_children=m_left->m_children+1;
369       m_height=m_left->m_height+1;
370       SIZELIMIT
371  }  }
372    
373    
# Line 385  DataLazy::DataLazy(DataAbstract_ptr left Line 377  DataLazy::DataLazy(DataAbstract_ptr left
377      m_op(op),      m_op(op),
378      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
379  {  {
380    cout << "Forming operator with " << left.get() << " " << right.get() << endl;
381     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
382     {     {
383      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
# Line 401  DataLazy::DataLazy(DataAbstract_ptr left Line 394  DataLazy::DataLazy(DataAbstract_ptr left
394     {     {
395      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
396      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
397    cout << "Right interpolation required " << right.get() << endl;
398     }     }
399     left->operandCheck(*right);     left->operandCheck(*right);
400    
401     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
402     {     {
403      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
404    cout << "Left is " << m_left->toString() << endl;
405     }     }
406     else     else
407     {     {
408      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
409    cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;
410     }     }
411     if (right->isLazy())     if (right->isLazy())
412     {     {
413      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
414    cout << "Right is " << m_right->toString() << endl;
415     }     }
416     else     else
417     {     {
418      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
419    cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;
420     }     }
421     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
422     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 437  DataLazy::DataLazy(DataAbstract_ptr left Line 435  DataLazy::DataLazy(DataAbstract_ptr left
435     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
436     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
437     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
438       m_children=m_left->m_children+m_right->m_children+2;
439       m_height=max(m_left->m_height,m_right->m_height)+1;
440       SIZELIMIT
441  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
442  }  }
443    
# Line 466  DataLazy::DataLazy(DataAbstract_ptr left Line 467  DataLazy::DataLazy(DataAbstract_ptr left
467      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
468      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
469     }     }
470     left->operandCheck(*right);  //    left->operandCheck(*right);
471    
472     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
473     {     {
# Line 501  DataLazy::DataLazy(DataAbstract_ptr left Line 502  DataLazy::DataLazy(DataAbstract_ptr left
502     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
503     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
504     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
505       m_children=m_left->m_children+m_right->m_children+2;
506       m_height=max(m_left->m_height,m_right->m_height)+1;
507       SIZELIMIT
508  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
509  }  }
510    
# Line 530  DataLazy::DataLazy(DataAbstract_ptr left Line 534  DataLazy::DataLazy(DataAbstract_ptr left
534     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
535     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
536     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
537       m_children=m_left->m_children+1;
538       m_height=m_left->m_height+1;
539       SIZELIMIT
540  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
541  }  }
542    
# Line 558  DataLazy::DataLazy(DataAbstract_ptr left Line 565  DataLazy::DataLazy(DataAbstract_ptr left
565     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
566     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
567     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
568       m_children=m_left->m_children+1;
569       m_height=m_left->m_height+1;
570       SIZELIMIT
571  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
572  }  }
573    
# Line 912  DataLazy::resolveNP1OUT(ValueType& v, si Line 922  DataLazy::resolveNP1OUT(ValueType& v, si
922    }    }
923      // 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
924    size_t subroffset=roffset+m_samplesize;    size_t subroffset=roffset+m_samplesize;
925    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);  LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
926      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
927    roffset=offset;    roffset=offset;
928      size_t loop=0;
929      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
930      size_t step=getNoValues();
931    switch (m_op)    switch (m_op)
932    {    {
933      case SYM:      case SYM:
934      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
935        {
936            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
937            subroffset+=step;
938            offset+=step;
939        }
940      break;      break;
941      case NSYM:      case NSYM:
942      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
943        {
944            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
945            subroffset+=step;
946            offset+=step;
947        }
948      break;      break;
949      default:      default:
950      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 975  DataLazy::resolveNP1OUT_P(ValueType& v,
975      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.");
976    }    }
977      // 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
978    size_t subroffset=roffset+m_samplesize;    size_t subroffset;
979    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);    const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
980    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
981    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
982    roffset=offset;    roffset=offset;
983      size_t loop=0;
984      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
985      size_t outstep=getNoValues();
986      size_t instep=m_left->getNoValues();
987    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
988    switch (m_op)    switch (m_op)
989    {    {
990      case TRACE:      case TRACE:
991           DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
992        {
993    size_t zz=sampleNo*getNumDPPSample()+loop;
994    if (zz==5800)
995    {
996    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
997    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
998    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
999    LAZYDEBUG(cerr << subroffset << endl;)
1000    LAZYDEBUG(cerr << "output=" << offset << endl;)
1001    }
1002                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1003    if (zz==5800)
1004    {
1005    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1006    }
1007            subroffset+=instep;
1008            offset+=outstep;
1009        }
1010      break;      break;
1011      case TRANS:      case TRANS:
1012           DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
1013        {
1014                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1015            subroffset+=instep;
1016            offset+=outstep;
1017        }
1018      break;      break;
1019      default:      default:
1020      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 1024  DataLazy::resolveNP1OUT_P(ValueType& v,
1024    
1025    
1026  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1027      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int j=0;j<onumsteps;++j)\
1028      { \      {\
1029         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1030         lroffset+=leftStep; \        { \
1031         rroffset+=rightStep; \  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1032    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1033             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1034    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1035             lroffset+=leftstep; \
1036             rroffset+=rightstep; \
1037          }\
1038          lroffset+=oleftstep;\
1039          rroffset+=orightstep;\
1040      }      }
1041    
1042  /*  /*
# Line 1013  LAZYDEBUG(cout << "Resolve binary: " << Line 1075  LAZYDEBUG(cout << "Resolve binary: " <<
1075    }    }
1076    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1077    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1078    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  
1079    {    {
1080      if (!leftScalar && !rightScalar)      throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1081      }
1082      size_t leftsize=m_left->getNoValues();
1083      size_t rightsize=m_right->getNoValues();
1084      size_t chunksize=1;           // how many doubles will be processed in one go
1085      int leftstep=0;       // how far should the left offset advance after each step
1086      int rightstep=0;
1087      int numsteps=0;       // total number of steps for the inner loop
1088      int oleftstep=0;  // the o variables refer to the outer loop
1089      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1090      int onumsteps=1;
1091      
1092      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1093      bool RES=(rightExp && rightScalar);
1094      bool LS=(!leftExp && leftScalar); // left is a single scalar
1095      bool RS=(!rightExp && rightScalar);
1096      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1097      bool RN=(!rightExp && !rightScalar);
1098      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1099      bool REN=(rightExp && !rightScalar);
1100    
1101      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1102      {
1103        chunksize=m_left->getNumDPPSample()*leftsize;
1104        leftstep=0;
1105        rightstep=0;
1106        numsteps=1;
1107      }
1108      else if (LES || RES)
1109      {
1110        chunksize=1;
1111        if (LES)        // left is an expanded scalar
1112      {      {
1113         throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");          if (RS)
1114            {
1115               leftstep=1;
1116               rightstep=0;
1117               numsteps=m_left->getNumDPPSample();
1118            }
1119            else        // RN or REN
1120            {
1121               leftstep=0;
1122               oleftstep=1;
1123               rightstep=1;
1124               orightstep=(RN ? -(int)rightsize : 0);
1125               numsteps=rightsize;
1126               onumsteps=m_left->getNumDPPSample();
1127            }
1128      }      }
1129      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());      else        // right is an expanded scalar
1130      chunksize=1;    // for scalar      {
1131    }              if (LS)
1132    int leftStep=((leftExp && (!rightExp || rightScalar))? m_right->getNoValues() : 0);          {
1133    int rightStep=((rightExp && (!leftExp || leftScalar))? m_left->getNoValues() : 0);             rightstep=1;
1134    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0             leftstep=0;
1135               numsteps=m_right->getNumDPPSample();
1136            }
1137            else
1138            {
1139               rightstep=0;
1140               orightstep=1;
1141               leftstep=1;
1142               oleftstep=(LN ? -(int)leftsize : 0);
1143               numsteps=leftsize;
1144               onumsteps=m_right->getNumDPPSample();
1145            }
1146        }
1147      }
1148      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1149      {
1150        if (LEN)    // and Right will be a single value
1151        {
1152            chunksize=rightsize;
1153            leftstep=rightsize;
1154            rightstep=0;
1155            numsteps=m_left->getNumDPPSample();
1156            if (RS)
1157            {
1158               numsteps*=leftsize;
1159            }
1160        }
1161        else    // REN
1162        {
1163            chunksize=leftsize;
1164            rightstep=leftsize;
1165            leftstep=0;
1166            numsteps=m_right->getNumDPPSample();
1167            if (LS)
1168            {
1169               numsteps*=rightsize;
1170            }
1171        }
1172      }
1173    
1174      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1175      // Get the values of sub-expressions      // Get the values of sub-expressions
1176    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1177    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
1178      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1179      // the right child starts further along.      // the right child starts further along.
1180    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1181    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1182    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1183    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1184    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1185    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1186    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1187    // LAZYDEBUG(
1188    // cout << "Results of bin" << endl;
1189    // cout << "Left=";
1190    // for (int i=lroffset;i<lroffset+m_left->m_samplesize;++i)
1191    // {
1192    // cout << (*left)[i] << " ";
1193    // }
1194    // cout << endl << "Right=";
1195    // for (int i=rroffset;i<rroffset+(m_right->m_readytype=='E'?m_right->m_samplesize:m_right->getNoValues());++i)
1196    // {
1197    // cout << (*right)[i] << " ";
1198    // }
1199    // cout << endl;
1200    // )
1201    
1202    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
1203    switch(m_op)    switch(m_op)
1204    {    {
# Line 1053  LAZYDEBUG(cout << "Resolve binary: " << Line 1220  LAZYDEBUG(cout << "Resolve binary: " <<
1220      default:      default:
1221      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1222    }    }
1223    roffset=offset;      roffset=offset;
1224    return &v;    return &v;
1225  }  }
1226    
1227    
1228    
1229  /*  /*
1230    \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.
1231    \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 1076  LAZYDEBUG(cout << "Resolve binary: " << Line 1244  LAZYDEBUG(cout << "Resolve binary: " <<
1244  DataTypes::ValueType*  DataTypes::ValueType*
1245  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
1246  {  {
1247  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1248    
1249    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
1250      // first work out which of the children are expanded      // first work out which of the children are expanded
1251    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1252    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1253    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1254    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);  /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1255    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1256    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1257      int rightStep=(rightExp?m_right->getNoValues() : 0);
1258    
1259      int resultStep=getNoValues();
1260    //   int resultStep=max(leftStep,rightStep);    // only one (at most) should be !=0
1261      // 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).
1262    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
1263    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);  
1264    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1265    
1266      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1267      gap+=m_right->getMaxSampleSize();
1268    
1269    
1270    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1271    
1272    
1273      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1274    
1275    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1276    cout << getNoValues() << endl;)
1277    LAZYDEBUG(cerr << "Result of left=";)
1278    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1279    for (int i=lroffset;i<lroffset+m_left->getNoValues();++i)
1280    {
1281    cout << (*left)[i] << " ";
1282    })
1283    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1284    LAZYDEBUG(cerr << "[" << rroffset << " .. " << rroffset+m_right->m_samplesize << "]" << endl;
1285    for (int i=rroffset;i<rroffset+m_right->m_samplesize;++i)
1286    {
1287    cerr << (*right)[i] << " ";
1288    }
1289    cerr << endl;
1290    )
1291    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1292    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1293    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1294    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1295    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1296    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1297    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1298    
1299    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
1300    switch(m_op)    switch(m_op)
1301    {    {
1302      case PROD:      case PROD:
1303      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1304      {      {
1305    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1306    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1307    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1308            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1309            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1310    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1311    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1312            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);
1313    LAZYDEBUG(cout << "Results--\n";
1314    for (int z=0;z<getNoValues();++z)
1315    {
1316    cout << resultp[z] << " ";
1317    }
1318    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1319    )
1320        lroffset+=leftStep;        lroffset+=leftStep;
1321        rroffset+=rightStep;        rroffset+=rightStep;
1322      }      }
# Line 1142  LAZYDEBUG(cout << "Resolve sample " << t Line 1361  LAZYDEBUG(cout << "Resolve sample " << t
1361      if (m_readytype=='C')      if (m_readytype=='C')
1362      {      {
1363      roffset=0;      roffset=0;
1364    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1365      return &(vec);      return &(vec);
1366      }      }
1367      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1368    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1369      return &(vec);      return &(vec);
1370    }    }
1371    if (m_readytype!='E')    if (m_readytype!='E')
# Line 1162  LAZYDEBUG(cout << "Resolve sample " << t Line 1383  LAZYDEBUG(cout << "Resolve sample " << t
1383    default:    default:
1384      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)+".");
1385    }    }
1386    
1387    }
1388    
1389    // This needs to do the work of the idenity constructor
1390    void
1391    DataLazy::resolveToIdentity()
1392    {
1393       if (m_op==IDENTITY)
1394        return;
1395       DataReady_ptr p=resolve();
1396       makeIdentity(p);
1397  }  }
1398    
1399    void DataLazy::makeIdentity(const DataReady_ptr& p)
1400    {
1401       m_op=IDENTITY;
1402       m_axis_offset=0;
1403       m_transpose=0;
1404       m_SL=m_SM=m_SR=0;
1405       m_children=m_height=0;
1406       m_id=p;
1407       if(p->isConstant()) {m_readytype='C';}
1408       else if(p->isExpanded()) {m_readytype='E';}
1409       else if (p->isTagged()) {m_readytype='T';}
1410       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1411       m_buffsRequired=1;
1412       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1413       m_maxsamplesize=m_samplesize;
1414       m_left.reset();
1415       m_right.reset();
1416    }
1417    
1418  // 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.
1419  // 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 1449  LAZYDEBUG(cout << "Buffer created with s
1449    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1450    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1451    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1452    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1453    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1454    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1455    {    {
1456          if (sample==0)  {ENABLEDEBUG}
1457    
1458    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1459  LAZYDEBUG(cout << "################################# " << sample << endl;)  LAZYDEBUG(cout << "################################# " << sample << endl;)
1460  #ifdef _OPENMP  #ifdef _OPENMP
1461      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 1463  LAZYDEBUG(cout << "#####################
1463      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.
1464  #endif  #endif
1465  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1466    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1467      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1468  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1469      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
1470      {      {
1471    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1472      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1473      }      }
1474    LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1475  LAZYDEBUG(cerr << "*********************************" << endl;)  LAZYDEBUG(cerr << "*********************************" << endl;)
1476        DISABLEDEBUG
1477    }    }
1478    return resptr;    return resptr;
1479  }  }
# Line 1233  DataLazy::toString() const Line 1491  DataLazy::toString() const
1491  void  void
1492  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1493  {  {
1494    //   oss << "[" << m_children <<";"<<m_height <<"]";
1495    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1496    {    {
1497    case G_IDENTITY:    case G_IDENTITY:

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

  ViewVC Help
Powered by ViewVC 1.1.26