/[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 2472 by jfenwick, Thu Jun 11 23:33:47 2009 UTC
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31    #include "EscriptParams.h"
32    
33    #include <iomanip>      // for some fancy formatting in debug
34    
35    // #define LAZYDEBUG(X) if (privdebug){X;}
36    #define LAZYDEBUG(X)
37    namespace
38    {
39    bool privdebug=false;
40    
41    #define ENABLEDEBUG privdebug=true;
42    #define DISABLEDEBUG privdebug=false;
43    }
44    
45    // #define SIZELIMIT
46    // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
47    //#define SIZELIMIT if ((m_height>7) || (m_children>127)) {resolveToIdentity();}
48    
49    #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}
50    
51  /*  /*
52  How does DataLazy work?  How does DataLazy work?
53  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 70  The convention that I use, is that the r Line 90  The convention that I use, is that the r
90  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.
91  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.
92    
93  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):
94  1) Add to the ES_optype.  1) Add to the ES_optype.
95  2) determine what opgroup your operation belongs to (X)  2) determine what opgroup your operation belongs to (X)
96  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 116  enum ES_opgroup
116     G_IDENTITY,     G_IDENTITY,
117     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
118     G_UNARY,     // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
119       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
120     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
121       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
122     G_TENSORPROD     // general tensor product     G_TENSORPROD     // general tensor product
123  };  };
124    
# Line 108  string ES_opstrings[]={"UNKNOWN","IDENTI Line 130  string ES_opstrings[]={"UNKNOWN","IDENTI
130              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
131              "asinh","acosh","atanh",              "asinh","acosh","atanh",
132              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
133              "1/","where>0","where<0","where>=0","where<=0",              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
134              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
135              "prod"};              "prod",
136  int ES_opcount=36;              "transpose", "trace"};
137    int ES_opcount=40;
138  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,
139              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
140              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
141              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
142              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
143              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
144              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
145              G_TENSORPROD};              G_TENSORPROD,
146                G_NP1OUT_P, G_NP1OUT_P};
147  inline  inline
148  ES_opgroup  ES_opgroup
149  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 177  resultShape(DataAbstract_ptr left, DataA Line 201  resultShape(DataAbstract_ptr left, DataA
201      return left->getShape();      return left->getShape();
202  }  }
203    
204    // return the shape for "op left"
205    
206    DataTypes::ShapeType
207    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
208    {
209        switch(op)
210        {
211            case TRANS:
212           {            // for the scoping of variables
213            const DataTypes::ShapeType& s=left->getShape();
214            DataTypes::ShapeType sh;
215            int rank=left->getRank();
216            if (axis_offset<0 || axis_offset>rank)
217            {
218                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
219                }
220                for (int i=0; i<rank; i++)
221            {
222               int index = (axis_offset+i)%rank;
223                   sh.push_back(s[index]); // Append to new shape
224                }
225            return sh;
226           }
227        break;
228        case TRACE:
229           {
230            int rank=left->getRank();
231            if (rank<2)
232            {
233               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
234            }
235            if ((axis_offset>rank-2) || (axis_offset<0))
236            {
237               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
238            }
239            if (rank==2)
240            {
241               return DataTypes::scalarShape;
242            }
243            else if (rank==3)
244            {
245               DataTypes::ShapeType sh;
246                   if (axis_offset==0)
247               {
248                    sh.push_back(left->getShape()[2]);
249                   }
250                   else     // offset==1
251               {
252                sh.push_back(left->getShape()[0]);
253                   }
254               return sh;
255            }
256            else if (rank==4)
257            {
258               DataTypes::ShapeType sh;
259               const DataTypes::ShapeType& s=left->getShape();
260                   if (axis_offset==0)
261               {
262                    sh.push_back(s[2]);
263                    sh.push_back(s[3]);
264                   }
265                   else if (axis_offset==1)
266               {
267                    sh.push_back(s[0]);
268                    sh.push_back(s[3]);
269                   }
270               else     // offset==2
271               {
272                sh.push_back(s[0]);
273                sh.push_back(s[1]);
274               }
275               return sh;
276            }
277            else        // unknown rank
278            {
279               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
280            }
281           }
282        break;
283            default:
284        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
285        }
286    }
287    
288  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
289  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
290  // 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 305  GTPShape(DataAbstract_ptr left, DataAbst
305    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)  { start1 = rank1-axis_offset; }
306    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"); }
307    
308      if (rank0<axis_offset)
309      {
310        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
311      }
312    
313    // Adjust the shapes for transpose    // Adjust the shapes for transpose
314    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 338  GTPShape(DataAbstract_ptr left, DataAbst
338       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
339       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
340    }    }
   return shape2;  
 }  
341    
342      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
343      {
344         ostringstream os;
345         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
346         throw DataException(os.str());
347      }
348    
349  // determine the number of points in the result of "left op right"    return shape2;
350  // 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)+".");  
 //    }  
 // }  
351    
352  // 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
353  // 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.
354  // The same goes for G_TENSORPROD  // The same goes for G_TENSORPROD
355    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
356    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
357    // multiple times
358  int  int
359  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
360  {  {
361     switch(getOpgroup(op))     switch(getOpgroup(op))
362     {     {
363     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
364     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
365     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY:
366       case G_UNARY_P: return max(left->getBuffsRequired(),1);
367     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
368       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
369     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
370     default:     default:
371      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
# Line 281  opToString(ES_optype op) Line 390  opToString(ES_optype op)
390    
391    
392  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
393      : 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)  
394  {  {
395     if (p->isLazy())     if (p->isLazy())
396     {     {
# Line 296  DataLazy::DataLazy(DataAbstract_ptr p) Line 401  DataLazy::DataLazy(DataAbstract_ptr p)
401     }     }
402     else     else
403     {     {
404      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
405      if(p->isConstant()) {m_readytype='C';}      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
406      else if(p->isExpanded()) {m_readytype='E';}      makeIdentity(dr);
407      else if (p->isTagged()) {m_readytype='T';}  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}  
408     }     }
409     m_buffsRequired=1;  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_maxsamplesize=m_samplesize;  
 cout << "(1)Lazy created with " << m_samplesize << endl;  
410  }  }
411    
412    
# Line 337  DataLazy::DataLazy(DataAbstract_ptr left Line 438  DataLazy::DataLazy(DataAbstract_ptr left
438     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
439     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
440     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
441       m_children=m_left->m_children+1;
442       m_height=m_left->m_height+1;
443       SIZELIMIT
444  }  }
445    
446    
# Line 346  DataLazy::DataLazy(DataAbstract_ptr left Line 450  DataLazy::DataLazy(DataAbstract_ptr left
450      m_op(op),      m_op(op),
451      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
452  {  {
453    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
454     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
455     {     {
456      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 362  DataLazy::DataLazy(DataAbstract_ptr left Line 467  DataLazy::DataLazy(DataAbstract_ptr left
467     {     {
468      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
469      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
470    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
471     }     }
472     left->operandCheck(*right);     left->operandCheck(*right);
473    
474     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
475     {     {
476      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
477    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
478     }     }
479     else     else
480     {     {
481      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
482    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
483     }     }
484     if (right->isLazy())     if (right->isLazy())
485     {     {
486      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
487    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
488     }     }
489     else     else
490     {     {
491      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
492    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
493     }     }
494     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
495     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 398  DataLazy::DataLazy(DataAbstract_ptr left Line 508  DataLazy::DataLazy(DataAbstract_ptr left
508     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
509     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
510     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
511  cout << "(3)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
512       m_height=max(m_left->m_height,m_right->m_height)+1;
513       SIZELIMIT
514    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
515  }  }
516    
517  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 427  DataLazy::DataLazy(DataAbstract_ptr left Line 540  DataLazy::DataLazy(DataAbstract_ptr left
540      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
541      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
542     }     }
543     left->operandCheck(*right);  //    left->operandCheck(*right);
544    
545     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
546     {     {
# Line 462  DataLazy::DataLazy(DataAbstract_ptr left Line 575  DataLazy::DataLazy(DataAbstract_ptr left
575     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
576     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
577     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
578  cout << "(4)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
579       m_height=max(m_left->m_height,m_right->m_height)+1;
580       SIZELIMIT
581    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
582    }
583    
584    
585    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
586        : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
587        m_op(op),
588        m_axis_offset(axis_offset),
589        m_transpose(0),
590        m_tol(0)
591    {
592       if ((getOpgroup(op)!=G_NP1OUT_P))
593       {
594        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
595       }
596       DataLazy_ptr lleft;
597       if (!left->isLazy())
598       {
599        lleft=DataLazy_ptr(new DataLazy(left));
600       }
601       else
602       {
603        lleft=dynamic_pointer_cast<DataLazy>(left);
604       }
605       m_readytype=lleft->m_readytype;
606       m_left=lleft;
607       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
608       m_samplesize=getNumDPPSample()*getNoValues();
609       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
610       m_children=m_left->m_children+1;
611       m_height=m_left->m_height+1;
612       SIZELIMIT
613    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
614  }  }
615    
616    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
617        : parent(left->getFunctionSpace(), left->getShape()),
618        m_op(op),
619        m_axis_offset(0),
620        m_transpose(0),
621        m_tol(tol)
622    {
623       if ((getOpgroup(op)!=G_UNARY_P))
624       {
625        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
626       }
627       DataLazy_ptr lleft;
628       if (!left->isLazy())
629       {
630        lleft=DataLazy_ptr(new DataLazy(left));
631       }
632       else
633       {
634        lleft=dynamic_pointer_cast<DataLazy>(left);
635       }
636       m_readytype=lleft->m_readytype;
637       m_left=lleft;
638       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
639       m_samplesize=getNumDPPSample()*getNoValues();
640       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
641       m_children=m_left->m_children+1;
642       m_height=m_left->m_height+1;
643       SIZELIMIT
644    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
645    }
646    
647  DataLazy::~DataLazy()  DataLazy::~DataLazy()
648  {  {
# Line 484  DataLazy::getMaxSampleSize() const Line 662  DataLazy::getMaxSampleSize() const
662      return m_maxsamplesize;      return m_maxsamplesize;
663  }  }
664    
665    
666    
667    size_t
668    DataLazy::getSampleBufferSize() const
669    {
670        return m_maxsamplesize*(max(1,m_buffsRequired));
671    }
672    
673  /*  /*
674    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
675    This does the work for the collapse method.    This does the work for the collapse method.
# Line 602  DataLazy::collapseToReady() Line 788  DataLazy::collapseToReady()
788      case LEZ:      case LEZ:
789      result=left.whereNonPositive();      result=left.whereNonPositive();
790      break;      break;
791        case NEZ:
792        result=left.whereNonZero(m_tol);
793        break;
794        case EZ:
795        result=left.whereZero(m_tol);
796        break;
797      case SYM:      case SYM:
798      result=left.symmetric();      result=left.symmetric();
799      break;      break;
# Line 611  DataLazy::collapseToReady() Line 803  DataLazy::collapseToReady()
803      case PROD:      case PROD:
804      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
805      break;      break;
806        case TRANS:
807        result=left.transpose(m_axis_offset);
808        break;
809        case TRACE:
810        result=left.trace(m_axis_offset);
811        break;
812      default:      default:
813      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)+".");
814    }    }
# Line 762  DataLazy::resolveUnary(ValueType& v, siz Line 960  DataLazy::resolveUnary(ValueType& v, siz
960      case LEZ:      case LEZ:
961      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));
962      break;      break;
963    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
964        case NEZ:
965        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
966        break;
967        case EZ:
968        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
969        break;
970    
971      default:      default:
972      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 975  DataLazy::resolveUnary(ValueType& v, siz
975  }  }
976    
977    
978    
979    
980    
981    
982  /*  /*
983    \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.
984    \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 794  DataLazy::resolveNP1OUT(ValueType& v, si Line 1003  DataLazy::resolveNP1OUT(ValueType& v, si
1003    }    }
1004      // 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
1005    size_t subroffset=roffset+m_samplesize;    size_t subroffset=roffset+m_samplesize;
1006    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);  LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
1007      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1008    roffset=offset;    roffset=offset;
1009      size_t loop=0;
1010      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1011      size_t step=getNoValues();
1012    switch (m_op)    switch (m_op)
1013    {    {
1014      case SYM:      case SYM:
1015      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1016        {
1017            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1018            subroffset+=step;
1019            offset+=step;
1020        }
1021      break;      break;
1022      case NSYM:      case NSYM:
1023      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1024        {
1025            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1026            subroffset+=step;
1027            offset+=step;
1028        }
1029      break;      break;
1030      default:      default:
1031      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 810  DataLazy::resolveNP1OUT(ValueType& v, si Line 1033  DataLazy::resolveNP1OUT(ValueType& v, si
1033    return &v;    return &v;
1034  }  }
1035    
1036    /*
1037      \brief Compute the value of the expression (unary operation) for the given sample.
1038      \return Vector which stores the value of the subexpression for the given sample.
1039      \param v A vector to store intermediate results.
1040      \param offset Index in v to begin storing results.
1041      \param sampleNo Sample number to evaluate.
1042      \param roffset (output parameter) the offset in the return vector where the result begins.
1043    
1044      The return value will be an existing vector so do not deallocate it.
1045      If the result is stored in v it should be stored at the offset given.
1046      Everything from offset to the end of v should be considered available for this method to use.
1047    */
1048    DataTypes::ValueType*
1049    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1050    {
1051        // we assume that any collapsing has been done before we get here
1052        // since we only have one argument we don't need to think about only
1053        // processing single points.
1054      if (m_readytype!='E')
1055      {
1056        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1057      }
1058        // since we can't write the result over the input, we need a result offset further along
1059      size_t subroffset;
1060      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1061    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1062    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1063      roffset=offset;
1064      size_t loop=0;
1065      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1066      size_t outstep=getNoValues();
1067      size_t instep=m_left->getNoValues();
1068    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1069      switch (m_op)
1070      {
1071        case TRACE:
1072        for (loop=0;loop<numsteps;++loop)
1073        {
1074    size_t zz=sampleNo*getNumDPPSample()+loop;
1075    if (zz==5800)
1076    {
1077    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
1078    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1079    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1080    LAZYDEBUG(cerr << subroffset << endl;)
1081    LAZYDEBUG(cerr << "output=" << offset << endl;)
1082    }
1083                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1084    if (zz==5800)
1085    {
1086    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1087    }
1088            subroffset+=instep;
1089            offset+=outstep;
1090        }
1091        break;
1092        case TRANS:
1093        for (loop=0;loop<numsteps;++loop)
1094        {
1095                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1096            subroffset+=instep;
1097            offset+=outstep;
1098        }
1099        break;
1100        default:
1101        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1102      }
1103      return &v;
1104    }
1105    
1106    
1107  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1108      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int j=0;j<onumsteps;++j)\
1109      { \      {\
1110         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1111         lroffset+=leftStep; \        { \
1112         rroffset+=rightStep; \  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1113    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1114             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1115    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1116             lroffset+=leftstep; \
1117             rroffset+=rightstep; \
1118          }\
1119          lroffset+=oleftstep;\
1120          rroffset+=orightstep;\
1121      }      }
1122    
1123  /*  /*
# Line 845  DataLazy::resolveNP1OUT(ValueType& v, si Line 1144  DataLazy::resolveNP1OUT(ValueType& v, si
1144  DataTypes::ValueType*  DataTypes::ValueType*
1145  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
1146  {  {
1147  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1148    
1149    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
1150      // first work out which of the children are expanded      // first work out which of the children are expanded
1151    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1152    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1153    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1154    int steps=(bigloops?1:getNumDPPSample());    {
1155    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'.");
1156    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1157    {    bool leftScalar=(m_left->getRank()==0);
1158      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1159      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1160      chunksize=1;    // for scalar    {
1161    }          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1162    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1163    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1164    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1165      size_t chunksize=1;           // how many doubles will be processed in one go
1166      int leftstep=0;       // how far should the left offset advance after each step
1167      int rightstep=0;
1168      int numsteps=0;       // total number of steps for the inner loop
1169      int oleftstep=0;  // the o variables refer to the outer loop
1170      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1171      int onumsteps=1;
1172      
1173      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1174      bool RES=(rightExp && rightScalar);
1175      bool LS=(!leftExp && leftScalar); // left is a single scalar
1176      bool RS=(!rightExp && rightScalar);
1177      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1178      bool RN=(!rightExp && !rightScalar);
1179      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1180      bool REN=(rightExp && !rightScalar);
1181    
1182      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1183      {
1184        chunksize=m_left->getNumDPPSample()*leftsize;
1185        leftstep=0;
1186        rightstep=0;
1187        numsteps=1;
1188      }
1189      else if (LES || RES)
1190      {
1191        chunksize=1;
1192        if (LES)        // left is an expanded scalar
1193        {
1194            if (RS)
1195            {
1196               leftstep=1;
1197               rightstep=0;
1198               numsteps=m_left->getNumDPPSample();
1199            }
1200            else        // RN or REN
1201            {
1202               leftstep=0;
1203               oleftstep=1;
1204               rightstep=1;
1205               orightstep=(RN ? -(int)rightsize : 0);
1206               numsteps=rightsize;
1207               onumsteps=m_left->getNumDPPSample();
1208            }
1209        }
1210        else        // right is an expanded scalar
1211        {
1212            if (LS)
1213            {
1214               rightstep=1;
1215               leftstep=0;
1216               numsteps=m_right->getNumDPPSample();
1217            }
1218            else
1219            {
1220               rightstep=0;
1221               orightstep=1;
1222               leftstep=1;
1223               oleftstep=(LN ? -(int)leftsize : 0);
1224               numsteps=leftsize;
1225               onumsteps=m_right->getNumDPPSample();
1226            }
1227        }
1228      }
1229      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1230      {
1231        if (LEN)    // and Right will be a single value
1232        {
1233            chunksize=rightsize;
1234            leftstep=rightsize;
1235            rightstep=0;
1236            numsteps=m_left->getNumDPPSample();
1237            if (RS)
1238            {
1239               numsteps*=leftsize;
1240            }
1241        }
1242        else    // REN
1243        {
1244            chunksize=leftsize;
1245            rightstep=leftsize;
1246            leftstep=0;
1247            numsteps=m_right->getNumDPPSample();
1248            if (LS)
1249            {
1250               numsteps*=rightsize;
1251            }
1252        }
1253      }
1254    
1255      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1256      // Get the values of sub-expressions      // Get the values of sub-expressions
1257    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1258    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
1259      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1260      // the right child starts further along.      // the right child starts further along.
1261    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1262    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1263    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1264    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1265    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1266    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1267    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1268    
1269    
1270    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
1271    switch(m_op)    switch(m_op)
1272    {    {
# Line 888  cout << "Resolve binary: " << toString() Line 1288  cout << "Resolve binary: " << toString()
1288      default:      default:
1289      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1290    }    }
1291    roffset=offset;      roffset=offset;
1292    return &v;    return &v;
1293  }  }
1294    
1295    
1296    
1297  /*  /*
1298    \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.
1299    \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 1312  cout << "Resolve binary: " << toString()
1312  DataTypes::ValueType*  DataTypes::ValueType*
1313  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
1314  {  {
1315  cout << "Resolve TensorProduct: " << toString() << endl;  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1316    
1317    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
1318      // first work out which of the children are expanded      // first work out which of the children are expanded
1319    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1320    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1321    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1322    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);  /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1323    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1324    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
1325      int rightStep=(rightExp?m_right->getNoValues() : 0);
1326    
1327      int resultStep=getNoValues();
1328      // 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).
1329    const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);    int gap=offset+m_samplesize;  
1330    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);  
1331    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1332    
1333      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1334      gap+=m_left->getMaxSampleSize();
1335    
1336    
1337    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1338    
1339    
1340      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1341    
1342    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1343    cout << getNoValues() << endl;)
1344    LAZYDEBUG(cerr << "Result of left=";)
1345    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1346    
1347    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1348    {
1349    cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1350    if (i%4==0) cout << endl;
1351    })
1352    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1353    LAZYDEBUG(
1354    for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1355    {
1356    cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1357    if (i%4==0) cout << endl;
1358    }
1359    cerr << endl;
1360    )
1361    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1362    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1363    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1364    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1365    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1366    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1367    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1368    
1369    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
1370    switch(m_op)    switch(m_op)
1371    {    {
1372      case PROD:      case PROD:
1373      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1374      {      {
1375    
1376    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1377    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1378    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1379    
1380            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1381            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1382    
1383    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1384    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1385    
1386            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);
1387    
1388    LAZYDEBUG(cout << "Results--\n";
1389    {
1390      DataVector dv(getNoValues());
1391    for (int z=0;z<getNoValues();++z)
1392    {
1393      cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1394      if (z%4==0) cout << endl;
1395      dv[z]=resultp[z];
1396    }
1397    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1398    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1399    }
1400    )
1401        lroffset+=leftStep;        lroffset+=leftStep;
1402        rroffset+=rightStep;        rroffset+=rightStep;
1403      }      }
# Line 965  cout << "Resolve TensorProduct: " << toS Line 1430  cout << "Resolve TensorProduct: " << toS
1430  const DataTypes::ValueType*  const DataTypes::ValueType*
1431  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1432  {  {
1433  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1434      // 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
1435    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1436    {    {
# Line 973  cout << "Resolve sample " << toString() Line 1438  cout << "Resolve sample " << toString()
1438    }    }
1439    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
1440    {    {
1441      const ValueType& vec=m_id->getVector();      const ValueType& vec=m_id->getVectorRO();
1442      if (m_readytype=='C')      if (m_readytype=='C')
1443      {      {
1444      roffset=0;      roffset=0;
1445    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1446      return &(vec);      return &(vec);
1447      }      }
1448      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1449    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1450      return &(vec);      return &(vec);
1451    }    }
1452    if (m_readytype!='E')    if (m_readytype!='E')
# Line 988  cout << "Resolve sample " << toString() Line 1455  cout << "Resolve sample " << toString()
1455    }    }
1456    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1457    {    {
1458    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY:
1459      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1460    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1461    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1462      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1463    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1464    default:    default:
1465      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)+".");
1466    }    }
1467    
1468    }
1469    
1470    const DataTypes::ValueType*
1471    DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
1472    {
1473    #ifdef _OPENMP
1474        int tid=omp_get_thread_num();
1475    #else
1476        int tid=0;
1477    #endif
1478        return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
1479  }  }
1480    
1481    
1482    // This needs to do the work of the idenity constructor
1483    void
1484    DataLazy::resolveToIdentity()
1485    {
1486       if (m_op==IDENTITY)
1487        return;
1488       DataReady_ptr p=resolve();
1489       makeIdentity(p);
1490    }
1491    
1492    void DataLazy::makeIdentity(const DataReady_ptr& p)
1493    {
1494       m_op=IDENTITY;
1495       m_axis_offset=0;
1496       m_transpose=0;
1497       m_SL=m_SM=m_SR=0;
1498       m_children=m_height=0;
1499       m_id=p;
1500       if(p->isConstant()) {m_readytype='C';}
1501       else if(p->isExpanded()) {m_readytype='E';}
1502       else if (p->isTagged()) {m_readytype='T';}
1503       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1504       m_buffsRequired=1;
1505       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1506       m_maxsamplesize=m_samplesize;
1507       m_left.reset();
1508       m_right.reset();
1509    }
1510    
1511  // 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.
1512  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
1513  DataReady_ptr  DataReady_ptr
1514  DataLazy::resolve()  DataLazy::resolve()
1515  {  {
1516    
1517  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1518  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
   
1519    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
1520    {    {
1521      collapse();      collapse();
# Line 1020  cout << "Buffers=" << m_buffsRequired << Line 1529  cout << "Buffers=" << m_buffsRequired <<
1529      // storage to evaluate its expression      // storage to evaluate its expression
1530    int numthreads=1;    int numthreads=1;
1531  #ifdef _OPENMP  #ifdef _OPENMP
1532    numthreads=getNumberOfThreads();    numthreads=omp_get_max_threads();
1533  #endif  #endif
1534    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1535  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1536    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1537    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
1538    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1539    int sample;    int sample;
1540    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
1541    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1542    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1543    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1544    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1545    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1546    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1547    {    {
1548  cout << "################################# " << sample << endl;        if (sample==0)  {ENABLEDEBUG}
1549    
1550    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1551    LAZYDEBUG(cout << "################################# " << sample << endl;)
1552  #ifdef _OPENMP  #ifdef _OPENMP
1553      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1554  #else  #else
1555      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.
1556  #endif  #endif
1557  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1558    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1559      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1560  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1561      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
1562      {      {
1563    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1564      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1565      }      }
1566  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1567    LAZYDEBUG(cerr << "*********************************" << endl;)
1568        DISABLEDEBUG
1569    }    }
1570    return resptr;    return resptr;
1571  }  }
# Line 1066  DataLazy::toString() const Line 1583  DataLazy::toString() const
1583  void  void
1584  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1585  {  {
1586    //    oss << "[" << m_children <<";"<<m_height <<"]";
1587    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1588    {    {
1589    case G_IDENTITY:    case G_IDENTITY:
# Line 1095  DataLazy::intoString(ostringstream& oss) Line 1613  DataLazy::intoString(ostringstream& oss)
1613      oss << ')';      oss << ')';
1614      break;      break;
1615    case G_UNARY:    case G_UNARY:
1616      case G_UNARY_P:
1617    case G_NP1OUT:    case G_NP1OUT:
1618      case G_NP1OUT_P:
1619      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1620      m_left->intoString(oss);      m_left->intoString(oss);
1621      oss << ')';      oss << ')';
# Line 1203  DataLazy::getPointOffset(int sampleNo, Line 1723  DataLazy::getPointOffset(int sampleNo,
1723    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");    throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1724  }  }
1725    
1726  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
1727  // to zero, all the tags from all the DataTags would be in the result.  // I have decided to let Data:: handle this issue.
 // However since they all have the same value (0) whether they are there or not should not matter.  
 // So I have decided that for all types this method will create a constant 0.  
 // It can be promoted up as required.  
 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management  
 // but we can deal with that if it arrises.  
1728  void  void
1729  DataLazy::setToZero()  DataLazy::setToZero()
1730  {  {
1731    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
1732    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1733    m_op=IDENTITY;  //   m_op=IDENTITY;
1734    m_right.reset();    //   m_right.reset();  
1735    m_left.reset();  //   m_left.reset();
1736    m_readytype='C';  //   m_readytype='C';
1737    m_buffsRequired=1;  //   m_buffsRequired=1;
1738    
1739      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1740    
1741    }
1742    
1743    bool
1744    DataLazy::actsExpanded() const
1745    {
1746        return (m_readytype=='E');
1747  }  }
1748    
1749  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26