/[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 2199 by jfenwick, Thu Jan 8 06:10: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    #include <iomanip>      // for some fancy formatting in debug
32    
33    // #define LAZYDEBUG(X) if (privdebug){X;}
34    #define LAZYDEBUG(X)
35    namespace
36    {
37    bool privdebug=false;
38    
39    #define ENABLEDEBUG privdebug=true;
40    #define DISABLEDEBUG privdebug=false;
41    }
42    
43    #define SIZELIMIT
44    // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
45    
46    
47  /*  /*
48  How does DataLazy work?  How does DataLazy work?
49  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 70  The convention that I use, is that the r Line 86  The convention that I use, is that the r
86  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.
87  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.
88    
89  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):
90  1) Add to the ES_optype.  1) Add to the ES_optype.
91  2) determine what opgroup your operation belongs to (X)  2) determine what opgroup your operation belongs to (X)
92  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 112  enum ES_opgroup
112     G_IDENTITY,     G_IDENTITY,
113     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
114     G_UNARY,     // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
115       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
116     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
117       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
118     G_TENSORPROD     // general tensor product     G_TENSORPROD     // general tensor product
119  };  };
120    
# Line 108  string ES_opstrings[]={"UNKNOWN","IDENTI Line 126  string ES_opstrings[]={"UNKNOWN","IDENTI
126              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
127              "asinh","acosh","atanh",              "asinh","acosh","atanh",
128              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
129              "1/","where>0","where<0","where>=0","where<=0",              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
130              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
131              "prod"};              "prod",
132  int ES_opcount=36;              "transpose", "trace"};
133    int ES_opcount=40;
134  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,
135              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
136              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
137              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
138              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
139              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
140              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
141              G_TENSORPROD};              G_TENSORPROD,
142                G_NP1OUT_P, G_NP1OUT_P};
143  inline  inline
144  ES_opgroup  ES_opgroup
145  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 177  resultShape(DataAbstract_ptr left, DataA Line 197  resultShape(DataAbstract_ptr left, DataA
197      return left->getShape();      return left->getShape();
198  }  }
199    
200    // return the shape for "op left"
201    
202    DataTypes::ShapeType
203    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
204    {
205        switch(op)
206        {
207            case TRANS:
208           {            // for the scoping of variables
209            const DataTypes::ShapeType& s=left->getShape();
210            DataTypes::ShapeType sh;
211            int rank=left->getRank();
212            if (axis_offset<0 || axis_offset>rank)
213            {
214                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
215                }
216                for (int i=0; i<rank; i++)
217            {
218               int index = (axis_offset+i)%rank;
219                   sh.push_back(s[index]); // Append to new shape
220                }
221            return sh;
222           }
223        break;
224        case TRACE:
225           {
226            int rank=left->getRank();
227            if (rank<2)
228            {
229               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
230            }
231            if ((axis_offset>rank-2) || (axis_offset<0))
232            {
233               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
234            }
235            if (rank==2)
236            {
237               return DataTypes::scalarShape;
238            }
239            else if (rank==3)
240            {
241               DataTypes::ShapeType sh;
242                   if (axis_offset==0)
243               {
244                    sh.push_back(left->getShape()[2]);
245                   }
246                   else     // offset==1
247               {
248                sh.push_back(left->getShape()[0]);
249                   }
250               return sh;
251            }
252            else if (rank==4)
253            {
254               DataTypes::ShapeType sh;
255               const DataTypes::ShapeType& s=left->getShape();
256                   if (axis_offset==0)
257               {
258                    sh.push_back(s[2]);
259                    sh.push_back(s[3]);
260                   }
261                   else if (axis_offset==1)
262               {
263                    sh.push_back(s[0]);
264                    sh.push_back(s[3]);
265                   }
266               else     // offset==2
267               {
268                sh.push_back(s[0]);
269                sh.push_back(s[1]);
270               }
271               return sh;
272            }
273            else        // unknown rank
274            {
275               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
276            }
277           }
278        break;
279            default:
280        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
281        }
282    }
283    
284  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
285  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
286  // 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 301  GTPShape(DataAbstract_ptr left, DataAbst
301    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)  { start1 = rank1-axis_offset; }
302    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"); }
303    
304      if (rank0<axis_offset)
305      {
306        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
307      }
308    
309    // Adjust the shapes for transpose    // Adjust the shapes for transpose
310    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 334  GTPShape(DataAbstract_ptr left, DataAbst
334       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
335       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
336    }    }
   return shape2;  
 }  
337    
338      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
339      {
340         ostringstream os;
341         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
342         throw DataException(os.str());
343      }
344    
345  // determine the number of points in the result of "left op right"    return shape2;
346  // 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)+".");  
 //    }  
 // }  
347    
348  // 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
349  // 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.
350  // The same goes for G_TENSORPROD  // The same goes for G_TENSORPROD
351    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefths.
352    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
353    // multiple times
354  int  int
355  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
356  {  {
357     switch(getOpgroup(op))     switch(getOpgroup(op))
358     {     {
359     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
360     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
361     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY:
362       case G_UNARY_P: return max(left->getBuffsRequired(),1);
363     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
364       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
365     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
366     default:     default:
367      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 386  opToString(ES_optype op)
386    
387    
388  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
389      : 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)  
390  {  {
391     if (p->isLazy())     if (p->isLazy())
392     {     {
# Line 296  DataLazy::DataLazy(DataAbstract_ptr p) Line 397  DataLazy::DataLazy(DataAbstract_ptr p)
397     }     }
398     else     else
399     {     {
400      m_id=dynamic_pointer_cast<DataReady>(p);      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
401      if(p->isConstant()) {m_readytype='C';}      makeIdentity(dr);
402      else if(p->isExpanded()) {m_readytype='E';}  LAZYDEBUG(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.");}  
403     }     }
404     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;  
405  }  }
406    
407    
# Line 337  DataLazy::DataLazy(DataAbstract_ptr left Line 433  DataLazy::DataLazy(DataAbstract_ptr left
433     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
434     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
435     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
436       m_children=m_left->m_children+1;
437       m_height=m_left->m_height+1;
438       SIZELIMIT
439  }  }
440    
441    
# Line 346  DataLazy::DataLazy(DataAbstract_ptr left Line 445  DataLazy::DataLazy(DataAbstract_ptr left
445      m_op(op),      m_op(op),
446      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
447  {  {
448    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
449     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
450     {     {
451      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 462  DataLazy::DataLazy(DataAbstract_ptr left
462     {     {
463      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
464      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
465    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
466     }     }
467     left->operandCheck(*right);     left->operandCheck(*right);
468    
469     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
470     {     {
471      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
472    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
473     }     }
474     else     else
475     {     {
476      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
477    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
478     }     }
479     if (right->isLazy())     if (right->isLazy())
480     {     {
481      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
482    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
483     }     }
484     else     else
485     {     {
486      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
487    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
488     }     }
489     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
490     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 398  DataLazy::DataLazy(DataAbstract_ptr left Line 503  DataLazy::DataLazy(DataAbstract_ptr left
503     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
504     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
505     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
506  cout << "(3)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
507       m_height=max(m_left->m_height,m_right->m_height)+1;
508       SIZELIMIT
509    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
510  }  }
511    
512  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 535  DataLazy::DataLazy(DataAbstract_ptr left
535      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
536      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
537     }     }
538     left->operandCheck(*right);  //    left->operandCheck(*right);
539    
540     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
541     {     {
# Line 462  DataLazy::DataLazy(DataAbstract_ptr left Line 570  DataLazy::DataLazy(DataAbstract_ptr left
570     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
571     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
572     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
573  cout << "(4)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
574       m_height=max(m_left->m_height,m_right->m_height)+1;
575       SIZELIMIT
576    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
577    }
578    
579    
580    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
581        : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
582        m_op(op),
583        m_axis_offset(axis_offset),
584        m_transpose(0),
585        m_tol(0)
586    {
587       if ((getOpgroup(op)!=G_NP1OUT_P))
588       {
589        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
590       }
591       DataLazy_ptr lleft;
592       if (!left->isLazy())
593       {
594        lleft=DataLazy_ptr(new DataLazy(left));
595       }
596       else
597       {
598        lleft=dynamic_pointer_cast<DataLazy>(left);
599       }
600       m_readytype=lleft->m_readytype;
601       m_left=lleft;
602       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
603       m_samplesize=getNumDPPSample()*getNoValues();
604       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
605       m_children=m_left->m_children+1;
606       m_height=m_left->m_height+1;
607       SIZELIMIT
608    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
609  }  }
610    
611    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
612        : parent(left->getFunctionSpace(), left->getShape()),
613        m_op(op),
614        m_axis_offset(0),
615        m_transpose(0),
616        m_tol(tol)
617    {
618       if ((getOpgroup(op)!=G_UNARY_P))
619       {
620        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
621       }
622       DataLazy_ptr lleft;
623       if (!left->isLazy())
624       {
625        lleft=DataLazy_ptr(new DataLazy(left));
626       }
627       else
628       {
629        lleft=dynamic_pointer_cast<DataLazy>(left);
630       }
631       m_readytype=lleft->m_readytype;
632       m_left=lleft;
633       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
634       m_samplesize=getNumDPPSample()*getNoValues();
635       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
636       m_children=m_left->m_children+1;
637       m_height=m_left->m_height+1;
638       SIZELIMIT
639    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
640    }
641    
642  DataLazy::~DataLazy()  DataLazy::~DataLazy()
643  {  {
# Line 602  DataLazy::collapseToReady() Line 775  DataLazy::collapseToReady()
775      case LEZ:      case LEZ:
776      result=left.whereNonPositive();      result=left.whereNonPositive();
777      break;      break;
778        case NEZ:
779        result=left.whereNonZero(m_tol);
780        break;
781        case EZ:
782        result=left.whereZero(m_tol);
783        break;
784      case SYM:      case SYM:
785      result=left.symmetric();      result=left.symmetric();
786      break;      break;
# Line 611  DataLazy::collapseToReady() Line 790  DataLazy::collapseToReady()
790      case PROD:      case PROD:
791      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
792      break;      break;
793        case TRANS:
794        result=left.transpose(m_axis_offset);
795        break;
796        case TRACE:
797        result=left.trace(m_axis_offset);
798        break;
799      default:      default:
800      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)+".");
801    }    }
# Line 762  DataLazy::resolveUnary(ValueType& v, siz Line 947  DataLazy::resolveUnary(ValueType& v, siz
947      case LEZ:      case LEZ:
948      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));
949      break;      break;
950    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
951        case NEZ:
952        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
953        break;
954        case EZ:
955        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
956        break;
957    
958      default:      default:
959      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 962  DataLazy::resolveUnary(ValueType& v, siz
962  }  }
963    
964    
965    
966    
967    
968    
969  /*  /*
970    \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.
971    \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 990  DataLazy::resolveNP1OUT(ValueType& v, si
990    }    }
991      // 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
992    size_t subroffset=roffset+m_samplesize;    size_t subroffset=roffset+m_samplesize;
993    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);  LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
994      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
995    roffset=offset;    roffset=offset;
996      size_t loop=0;
997      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
998      size_t step=getNoValues();
999    switch (m_op)    switch (m_op)
1000    {    {
1001      case SYM:      case SYM:
1002      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1003        {
1004            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1005            subroffset+=step;
1006            offset+=step;
1007        }
1008      break;      break;
1009      case NSYM:      case NSYM:
1010      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1011        {
1012            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1013            subroffset+=step;
1014            offset+=step;
1015        }
1016      break;      break;
1017      default:      default:
1018      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 1020  DataLazy::resolveNP1OUT(ValueType& v, si
1020    return &v;    return &v;
1021  }  }
1022    
1023    /*
1024      \brief Compute the value of the expression (unary operation) for the given sample.
1025      \return Vector which stores the value of the subexpression for the given sample.
1026      \param v A vector to store intermediate results.
1027      \param offset Index in v to begin storing results.
1028      \param sampleNo Sample number to evaluate.
1029      \param roffset (output parameter) the offset in the return vector where the result begins.
1030    
1031      The return value will be an existing vector so do not deallocate it.
1032      If the result is stored in v it should be stored at the offset given.
1033      Everything from offset to the end of v should be considered available for this method to use.
1034    */
1035    DataTypes::ValueType*
1036    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1037    {
1038        // we assume that any collapsing has been done before we get here
1039        // since we only have one argument we don't need to think about only
1040        // processing single points.
1041      if (m_readytype!='E')
1042      {
1043        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1044      }
1045        // since we can't write the result over the input, we need a result offset further along
1046      size_t subroffset;
1047      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1048    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1049    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1050      roffset=offset;
1051      size_t loop=0;
1052      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1053      size_t outstep=getNoValues();
1054      size_t instep=m_left->getNoValues();
1055    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1056      switch (m_op)
1057      {
1058        case TRACE:
1059        for (loop=0;loop<numsteps;++loop)
1060        {
1061    size_t zz=sampleNo*getNumDPPSample()+loop;
1062    if (zz==5800)
1063    {
1064    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
1065    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1066    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1067    LAZYDEBUG(cerr << subroffset << endl;)
1068    LAZYDEBUG(cerr << "output=" << offset << endl;)
1069    }
1070                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1071    if (zz==5800)
1072    {
1073    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1074    }
1075            subroffset+=instep;
1076            offset+=outstep;
1077        }
1078        break;
1079        case TRANS:
1080        for (loop=0;loop<numsteps;++loop)
1081        {
1082                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1083            subroffset+=instep;
1084            offset+=outstep;
1085        }
1086        break;
1087        default:
1088        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1089      }
1090      return &v;
1091    }
1092    
1093    
1094  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1095      for (int i=0;i<steps;++i,resultp+=resultStep) \      for (int j=0;j<onumsteps;++j)\
1096      { \      {\
1097         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1098         lroffset+=leftStep; \        { \
1099         rroffset+=rightStep; \  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1100    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1101             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1102    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1103             lroffset+=leftstep; \
1104             rroffset+=rightstep; \
1105          }\
1106          lroffset+=oleftstep;\
1107          rroffset+=orightstep;\
1108      }      }
1109    
1110  /*  /*
# Line 845  DataLazy::resolveNP1OUT(ValueType& v, si Line 1131  DataLazy::resolveNP1OUT(ValueType& v, si
1131  DataTypes::ValueType*  DataTypes::ValueType*
1132  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
1133  {  {
1134  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1135    
1136    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
1137      // first work out which of the children are expanded      // first work out which of the children are expanded
1138    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1139    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1140    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1141    int steps=(bigloops?1:getNumDPPSample());    {
1142    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'.");
1143    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1144    {    bool leftScalar=(m_left->getRank()==0);
1145      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1146      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1147      chunksize=1;    // for scalar    {
1148    }          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1149    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1150    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1151    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1152      size_t chunksize=1;           // how many doubles will be processed in one go
1153      int leftstep=0;       // how far should the left offset advance after each step
1154      int rightstep=0;
1155      int numsteps=0;       // total number of steps for the inner loop
1156      int oleftstep=0;  // the o variables refer to the outer loop
1157      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1158      int onumsteps=1;
1159      
1160      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1161      bool RES=(rightExp && rightScalar);
1162      bool LS=(!leftExp && leftScalar); // left is a single scalar
1163      bool RS=(!rightExp && rightScalar);
1164      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1165      bool RN=(!rightExp && !rightScalar);
1166      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1167      bool REN=(rightExp && !rightScalar);
1168    
1169      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1170      {
1171        chunksize=m_left->getNumDPPSample()*leftsize;
1172        leftstep=0;
1173        rightstep=0;
1174        numsteps=1;
1175      }
1176      else if (LES || RES)
1177      {
1178        chunksize=1;
1179        if (LES)        // left is an expanded scalar
1180        {
1181            if (RS)
1182            {
1183               leftstep=1;
1184               rightstep=0;
1185               numsteps=m_left->getNumDPPSample();
1186            }
1187            else        // RN or REN
1188            {
1189               leftstep=0;
1190               oleftstep=1;
1191               rightstep=1;
1192               orightstep=(RN ? -(int)rightsize : 0);
1193               numsteps=rightsize;
1194               onumsteps=m_left->getNumDPPSample();
1195            }
1196        }
1197        else        // right is an expanded scalar
1198        {
1199            if (LS)
1200            {
1201               rightstep=1;
1202               leftstep=0;
1203               numsteps=m_right->getNumDPPSample();
1204            }
1205            else
1206            {
1207               rightstep=0;
1208               orightstep=1;
1209               leftstep=1;
1210               oleftstep=(LN ? -(int)leftsize : 0);
1211               numsteps=leftsize;
1212               onumsteps=m_right->getNumDPPSample();
1213            }
1214        }
1215      }
1216      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1217      {
1218        if (LEN)    // and Right will be a single value
1219        {
1220            chunksize=rightsize;
1221            leftstep=rightsize;
1222            rightstep=0;
1223            numsteps=m_left->getNumDPPSample();
1224            if (RS)
1225            {
1226               numsteps*=leftsize;
1227            }
1228        }
1229        else    // REN
1230        {
1231            chunksize=leftsize;
1232            rightstep=leftsize;
1233            leftstep=0;
1234            numsteps=m_right->getNumDPPSample();
1235            if (LS)
1236            {
1237               numsteps*=rightsize;
1238            }
1239        }
1240      }
1241    
1242      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1243      // Get the values of sub-expressions      // Get the values of sub-expressions
1244    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1245    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
1246      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1247      // the right child starts further along.      // the right child starts further along.
1248    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1249    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1250    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1251    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1252    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1253    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1254    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1255    
1256    
1257    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
1258    switch(m_op)    switch(m_op)
1259    {    {
# Line 888  cout << "Resolve binary: " << toString() Line 1275  cout << "Resolve binary: " << toString()
1275      default:      default:
1276      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1277    }    }
1278    roffset=offset;      roffset=offset;
1279    return &v;    return &v;
1280  }  }
1281    
1282    
1283    
1284  /*  /*
1285    \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.
1286    \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 1299  cout << "Resolve binary: " << toString()
1299  DataTypes::ValueType*  DataTypes::ValueType*
1300  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
1301  {  {
1302  cout << "Resolve TensorProduct: " << toString() << endl;  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1303    
1304    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
1305      // first work out which of the children are expanded      // first work out which of the children are expanded
1306    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1307    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1308    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1309    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);  /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1310    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1311    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
1312      int rightStep=(rightExp?m_right->getNoValues() : 0);
1313    
1314      int resultStep=getNoValues();
1315      // 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).
1316    const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);    int gap=offset+m_samplesize;  
1317    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);  
1318    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1319    
1320      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1321      gap+=m_left->getMaxSampleSize();
1322    
1323    
1324    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1325    
1326    
1327      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1328    
1329    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1330    cout << getNoValues() << endl;)
1331    LAZYDEBUG(cerr << "Result of left=";)
1332    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1333    
1334    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1335    {
1336    cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1337    if (i%4==0) cout << endl;
1338    })
1339    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1340    LAZYDEBUG(
1341    for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1342    {
1343    cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1344    if (i%4==0) cout << endl;
1345    }
1346    cerr << endl;
1347    )
1348    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1349    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1350    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1351    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1352    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1353    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1354    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1355    
1356    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
1357    switch(m_op)    switch(m_op)
1358    {    {
1359      case PROD:      case PROD:
1360      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1361      {      {
1362    
1363    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1364    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1365    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1366    
1367            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1368            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1369    
1370    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1371    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1372    
1373            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);
1374    
1375    LAZYDEBUG(cout << "Results--\n";
1376    {
1377      DataVector dv(getNoValues());
1378    for (int z=0;z<getNoValues();++z)
1379    {
1380      cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1381      if (z%4==0) cout << endl;
1382      dv[z]=resultp[z];
1383    }
1384    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1385    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1386    }
1387    )
1388        lroffset+=leftStep;        lroffset+=leftStep;
1389        rroffset+=rightStep;        rroffset+=rightStep;
1390      }      }
# Line 965  cout << "Resolve TensorProduct: " << toS Line 1417  cout << "Resolve TensorProduct: " << toS
1417  const DataTypes::ValueType*  const DataTypes::ValueType*
1418  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1419  {  {
1420  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1421      // 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
1422    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1423    {    {
# Line 977  cout << "Resolve sample " << toString() Line 1429  cout << "Resolve sample " << toString()
1429      if (m_readytype=='C')      if (m_readytype=='C')
1430      {      {
1431      roffset=0;      roffset=0;
1432    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1433      return &(vec);      return &(vec);
1434      }      }
1435      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1436    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1437      return &(vec);      return &(vec);
1438    }    }
1439    if (m_readytype!='E')    if (m_readytype!='E')
# Line 988  cout << "Resolve sample " << toString() Line 1442  cout << "Resolve sample " << toString()
1442    }    }
1443    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1444    {    {
1445    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY:
1446      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1447    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1448    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1449      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1450    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1451    default:    default:
1452      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)+".");
1453    }    }
1454    
1455  }  }
1456    
1457    // This needs to do the work of the idenity constructor
1458    void
1459    DataLazy::resolveToIdentity()
1460    {
1461       if (m_op==IDENTITY)
1462        return;
1463       DataReady_ptr p=resolve();
1464       makeIdentity(p);
1465    }
1466    
1467    void DataLazy::makeIdentity(const DataReady_ptr& p)
1468    {
1469       m_op=IDENTITY;
1470       m_axis_offset=0;
1471       m_transpose=0;
1472       m_SL=m_SM=m_SR=0;
1473       m_children=m_height=0;
1474       m_id=p;
1475       if(p->isConstant()) {m_readytype='C';}
1476       else if(p->isExpanded()) {m_readytype='E';}
1477       else if (p->isTagged()) {m_readytype='T';}
1478       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1479       m_buffsRequired=1;
1480       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1481       m_maxsamplesize=m_samplesize;
1482       m_left.reset();
1483       m_right.reset();
1484    }
1485    
1486  // 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.
1487  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
# Line 1004  DataReady_ptr Line 1489  DataReady_ptr
1489  DataLazy::resolve()  DataLazy::resolve()
1490  {  {
1491    
1492  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1493  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1494    
1495    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
1496    {    {
# Line 1023  cout << "Buffers=" << m_buffsRequired << Line 1508  cout << "Buffers=" << m_buffsRequired <<
1508    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
1509  #endif  #endif
1510    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1511  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1512    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1513    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1514    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
# Line 1032  cout << "Buffer created with size=" << v Line 1517  cout << "Buffer created with size=" << v
1517    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1518    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Vector storing the answer
1519    size_t resoffset=0;       // where in the vector to find the answer    size_t resoffset=0;       // where in the vector to find the answer
1520    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1521    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1522    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1523    {    {
1524  cout << "################################# " << sample << endl;        if (sample==0)  {ENABLEDEBUG}
1525    
1526    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1527    LAZYDEBUG(cout << "################################# " << sample << endl;)
1528  #ifdef _OPENMP  #ifdef _OPENMP
1529      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1530  #else  #else
1531      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.
1532  #endif  #endif
1533  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1534    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1535      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1536  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1537      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
1538      {      {
1539    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1540      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1541      }      }
1542  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1543    LAZYDEBUG(cerr << "*********************************" << endl;)
1544        DISABLEDEBUG
1545    }    }
1546    return resptr;    return resptr;
1547  }  }
# Line 1066  DataLazy::toString() const Line 1559  DataLazy::toString() const
1559  void  void
1560  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1561  {  {
1562    //   oss << "[" << m_children <<";"<<m_height <<"]";
1563    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1564    {    {
1565    case G_IDENTITY:    case G_IDENTITY:
# Line 1095  DataLazy::intoString(ostringstream& oss) Line 1589  DataLazy::intoString(ostringstream& oss)
1589      oss << ')';      oss << ')';
1590      break;      break;
1591    case G_UNARY:    case G_UNARY:
1592      case G_UNARY_P:
1593    case G_NP1OUT:    case G_NP1OUT:
1594      case G_NP1OUT_P:
1595      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1596      m_left->intoString(oss);      m_left->intoString(oss);
1597      oss << ')';      oss << ')';

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

  ViewVC Help
Powered by ViewVC 1.1.26