/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/clazy/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2147 by jfenwick, Wed Dec 10 04:41:26 2008 UTC revision 2791 by jfenwick, Tue Dec 1 03:36:36 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31  // #define LAZYDEBUG(X) X;  #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)  #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 if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl;resolveToIdentity();}
46    
47    #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}
48    
49    
50  /*  /*
51  How does DataLazy work?  How does DataLazy work?
# Line 93  namespace escript Line 109  namespace escript
109  namespace  namespace
110  {  {
111    
112    
113    // enabling this will print out when ever the maximum stacksize used by resolve increases
114    // it assumes _OPENMP is also in use
115    //#define LAZY_STACK_PROF
116    
117    
118    
119    #ifndef _OPENMP
120      #ifdef LAZY_STACK_PROF
121      #undef LAZY_STACK_PROF
122      #endif
123    #endif
124    
125    
126    #ifdef LAZY_STACK_PROF
127    std::vector<void*> stackstart(getNumberOfThreads());
128    std::vector<void*> stackend(getNumberOfThreads());
129    size_t maxstackuse=0;
130    #endif
131    
132  enum ES_opgroup  enum ES_opgroup
133  {  {
134     G_UNKNOWN,     G_UNKNOWN,
# Line 102  enum ES_opgroup Line 138  enum ES_opgroup
138     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
139     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
140     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
141     G_TENSORPROD     // general tensor product     G_TENSORPROD,    // general tensor product
142       G_NP1OUT_2P,     // non-pointwise op with one output requiring two params
143       G_REDUCTION      // non-pointwise unary op with a scalar output
144  };  };
145    
146    
# Line 116  string ES_opstrings[]={"UNKNOWN","IDENTI Line 154  string ES_opstrings[]={"UNKNOWN","IDENTI
154              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
155              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
156              "prod",              "prod",
157              "transpose", "trace"};              "transpose", "trace",
158  int ES_opcount=40;              "swapaxes",
159                "minval", "maxval"};
160    int ES_opcount=43;
161  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,
162              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
163              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
# Line 126  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT Line 166  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT
166              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35
167              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
168              G_TENSORPROD,              G_TENSORPROD,
169              G_NP1OUT_P, G_NP1OUT_P};              G_NP1OUT_P, G_NP1OUT_P,
170                G_NP1OUT_2P,
171                G_REDUCTION, G_REDUCTION};
172  inline  inline
173  ES_opgroup  ES_opgroup
174  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 171  resultShape(DataAbstract_ptr left, DataA Line 213  resultShape(DataAbstract_ptr left, DataA
213        {        {
214          throw DataException("Shapes not the name - shapes must match for (point)binary operations.");          throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
215        }        }
216    
217        if (left->getRank()==0)   // we need to allow scalar * anything        if (left->getRank()==0)   // we need to allow scalar * anything
218        {        {
219          return right->getShape();          return right->getShape();
# Line 187  resultShape(DataAbstract_ptr left, DataA Line 230  resultShape(DataAbstract_ptr left, DataA
230  // return the shape for "op left"  // return the shape for "op left"
231    
232  DataTypes::ShapeType  DataTypes::ShapeType
233  resultShape(DataAbstract_ptr left, ES_optype op)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
234  {  {
235      switch(op)      switch(op)
236      {      {
237          case TRANS:          case TRANS:
238          return left->getShape();         {            // for the scoping of variables
239            const DataTypes::ShapeType& s=left->getShape();
240            DataTypes::ShapeType sh;
241            int rank=left->getRank();
242            if (axis_offset<0 || axis_offset>rank)
243            {
244                stringstream e;
245                e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
246                throw DataException(e.str());
247            }
248            for (int i=0; i<rank; i++)
249            {
250               int index = (axis_offset+i)%rank;
251               sh.push_back(s[index]); // Append to new shape
252            }
253            return sh;
254           }
255      break;      break;
256      case TRACE:      case TRACE:
257          return DataTypes::scalarShape;         {
258            int rank=left->getRank();
259            if (rank<2)
260            {
261               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
262            }
263            if ((axis_offset>rank-2) || (axis_offset<0))
264            {
265               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
266            }
267            if (rank==2)
268            {
269               return DataTypes::scalarShape;
270            }
271            else if (rank==3)
272            {
273               DataTypes::ShapeType sh;
274                   if (axis_offset==0)
275               {
276                    sh.push_back(left->getShape()[2]);
277                   }
278                   else     // offset==1
279               {
280                sh.push_back(left->getShape()[0]);
281                   }
282               return sh;
283            }
284            else if (rank==4)
285            {
286               DataTypes::ShapeType sh;
287               const DataTypes::ShapeType& s=left->getShape();
288                   if (axis_offset==0)
289               {
290                    sh.push_back(s[2]);
291                    sh.push_back(s[3]);
292                   }
293                   else if (axis_offset==1)
294               {
295                    sh.push_back(s[0]);
296                    sh.push_back(s[3]);
297                   }
298               else     // offset==2
299               {
300                sh.push_back(s[0]);
301                sh.push_back(s[1]);
302               }
303               return sh;
304            }
305            else        // unknown rank
306            {
307               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
308            }
309           }
310      break;      break;
311          default:          default:
312      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
313      }      }
314  }  }
315    
316    DataTypes::ShapeType
317    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
318    {
319         // This code taken from the Data.cpp swapaxes() method
320         // Some of the checks are probably redundant here
321         int axis0_tmp,axis1_tmp;
322         const DataTypes::ShapeType& s=left->getShape();
323         DataTypes::ShapeType out_shape;
324         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
325         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
326         int rank=left->getRank();
327         if (rank<2) {
328            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
329         }
330         if (axis0<0 || axis0>rank-1) {
331            stringstream e;
332            e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
333            throw DataException(e.str());
334         }
335         if (axis1<0 || axis1>rank-1) {
336            stringstream e;
337            e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
338            throw DataException(e.str());
339         }
340         if (axis0 == axis1) {
341             throw DataException("Error - Data::swapaxes: axis indices must be different.");
342         }
343         if (axis0 > axis1) {
344             axis0_tmp=axis1;
345             axis1_tmp=axis0;
346         } else {
347             axis0_tmp=axis0;
348             axis1_tmp=axis1;
349         }
350         for (int i=0; i<rank; i++) {
351           if (i == axis0_tmp) {
352              out_shape.push_back(s[axis1_tmp]);
353           } else if (i == axis1_tmp) {
354              out_shape.push_back(s[axis0_tmp]);
355           } else {
356              out_shape.push_back(s[i]);
357           }
358         }
359        return out_shape;
360    }
361    
362    
363  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
364  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
365  // 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 266  GTPShape(DataAbstract_ptr left, DataAbst Line 424  GTPShape(DataAbstract_ptr left, DataAbst
424    return shape2;    return shape2;
425  }  }
426    
   
 // determine the number of points in the result of "left op right"  
 // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here  
 // size_t  
 // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
 // {  
 //    switch (getOpgroup(op))  
 //    {  
 //    case G_BINARY: return left->getLength();  
 //    case G_UNARY: return left->getLength();  
 //    case G_NP1OUT: return left->getLength();  
 //    default:  
 //  throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");  
 //    }  
 // }  
   
 // determine the number of samples requires to evaluate an expression combining left and right  
 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.  
 // The same goes for G_TENSORPROD  
 int  
 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  
 {  
    switch(getOpgroup(op))  
    {  
    case G_IDENTITY: return 1;  
    case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);  
    case G_UNARY:  
    case G_UNARY_P: return max(left->getBuffsRequired(),1);  
    case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);  
    case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);  
    case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);  
    default:  
     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");  
    }  
 }  
   
   
427  }   // end anonymous namespace  }   // end anonymous namespace
428    
429    
# Line 318  opToString(ES_optype op) Line 439  opToString(ES_optype op)
439    return ES_opstrings[op];    return ES_opstrings[op];
440  }  }
441    
442    void DataLazy::LazyNodeSetup()
443    {
444    #ifdef _OPENMP
445        int numthreads=omp_get_max_threads();
446        m_samples.resize(numthreads*m_samplesize);
447        m_sampleids=new int[numthreads];
448        for (int i=0;i<numthreads;++i)
449        {
450            m_sampleids[i]=-1;  
451        }
452    #else
453        m_samples.resize(m_samplesize);
454        m_sampleids=new int[1];
455        m_sampleids[0]=-1;
456    #endif  // _OPENMP
457    }
458    
459    
460    // Creates an identity node
461  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
462      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
463      m_op(IDENTITY),      ,m_sampleids(0),
464      m_axis_offset(0),      m_samples(1)
     m_transpose(0),  
     m_SL(0), m_SM(0), m_SR(0)  
465  {  {
466     if (p->isLazy())     if (p->isLazy())
467     {     {
# Line 335  DataLazy::DataLazy(DataAbstract_ptr p) Line 472  DataLazy::DataLazy(DataAbstract_ptr p)
472     }     }
473     else     else
474     {     {
475      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
476      if(p->isConstant()) {m_readytype='C';}      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
477      else if(p->isExpanded()) {m_readytype='E';}      makeIdentity(dr);
478      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.");}  
479     }     }
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_maxsamplesize=m_samplesize;  
480  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
481  }  }
482    
   
   
   
483  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
484      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
485      m_op(op),      m_op(op),
486      m_axis_offset(0),      m_axis_offset(0),
487      m_transpose(0),      m_transpose(0),
488      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
489  {  {
490     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
491     {     {
492      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
493     }     }
# Line 373  DataLazy::DataLazy(DataAbstract_ptr left Line 503  DataLazy::DataLazy(DataAbstract_ptr left
503     }     }
504     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
505     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
506     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
507     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
508       m_height=m_left->m_height+1;
509       LazyNodeSetup();
510       SIZELIMIT
511  }  }
512    
513    
# Line 385  DataLazy::DataLazy(DataAbstract_ptr left Line 517  DataLazy::DataLazy(DataAbstract_ptr left
517      m_op(op),      m_op(op),
518      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
519  {  {
520    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
521     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
522     {     {
523      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
# Line 401  DataLazy::DataLazy(DataAbstract_ptr left Line 534  DataLazy::DataLazy(DataAbstract_ptr left
534     {     {
535      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
536      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
537    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
538     }     }
539     left->operandCheck(*right);     left->operandCheck(*right);
540    
541     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
542     {     {
543      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
544    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
545     }     }
546     else     else
547     {     {
548      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
549    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
550     }     }
551     if (right->isLazy())     if (right->isLazy())
552     {     {
553      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
554    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
555     }     }
556     else     else
557     {     {
558      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
559    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
560     }     }
561     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
562     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 435  DataLazy::DataLazy(DataAbstract_ptr left Line 573  DataLazy::DataLazy(DataAbstract_ptr left
573      m_readytype='C';      m_readytype='C';
574     }     }
575     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
576     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_children=m_left->m_children+m_right->m_children+2;
577     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
578       LazyNodeSetup();
579       SIZELIMIT
580  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
581  }  }
582    
# Line 466  DataLazy::DataLazy(DataAbstract_ptr left Line 606  DataLazy::DataLazy(DataAbstract_ptr left
606      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
607      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
608     }     }
609     left->operandCheck(*right);  //    left->operandCheck(*right);
610    
611     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
612     {     {
# Line 499  DataLazy::DataLazy(DataAbstract_ptr left Line 639  DataLazy::DataLazy(DataAbstract_ptr left
639      m_readytype='C';      m_readytype='C';
640     }     }
641     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
642     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_children=m_left->m_children+m_right->m_children+2;
643     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
644       LazyNodeSetup();
645       SIZELIMIT
646  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
647  }  }
648    
649    
650  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
651      : parent(left->getFunctionSpace(), resultShape(left,op)),      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
652      m_op(op),      m_op(op),
653      m_axis_offset(axis_offset),      m_axis_offset(axis_offset),
654      m_transpose(0),      m_transpose(0),
# Line 527  DataLazy::DataLazy(DataAbstract_ptr left Line 669  DataLazy::DataLazy(DataAbstract_ptr left
669     }     }
670     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
671     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
672     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
673     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
674       m_height=m_left->m_height+1;
675       LazyNodeSetup();
676       SIZELIMIT
677  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
678  }  }
679    
# Line 555  DataLazy::DataLazy(DataAbstract_ptr left Line 699  DataLazy::DataLazy(DataAbstract_ptr left
699     }     }
700     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
701     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
702     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
703     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
704       m_height=m_left->m_height+1;
705       LazyNodeSetup();
706       SIZELIMIT
707  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
708  }  }
709    
 DataLazy::~DataLazy()  
 {  
 }  
   
710    
711  int  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
712  DataLazy::getBuffsRequired() const      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
713        m_op(op),
714        m_axis_offset(axis0),
715        m_transpose(axis1),
716        m_tol(0)
717  {  {
718      return m_buffsRequired;     if ((getOpgroup(op)!=G_NP1OUT_2P))
719       {
720        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
721       }
722       DataLazy_ptr lleft;
723       if (!left->isLazy())
724       {
725        lleft=DataLazy_ptr(new DataLazy(left));
726       }
727       else
728       {
729        lleft=dynamic_pointer_cast<DataLazy>(left);
730       }
731       m_readytype=lleft->m_readytype;
732       m_left=lleft;
733       m_samplesize=getNumDPPSample()*getNoValues();
734       m_children=m_left->m_children+1;
735       m_height=m_left->m_height+1;
736       LazyNodeSetup();
737       SIZELIMIT
738    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
739  }  }
740    
741    DataLazy::~DataLazy()
 size_t  
 DataLazy::getMaxSampleSize() const  
742  {  {
743      return m_maxsamplesize;     delete[] m_sampleids;
744  }  }
745    
746    
747  /*  /*
748    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
749    This does the work for the collapse method.    This does the work for the collapse method.
# Line 718  DataLazy::collapseToReady() Line 883  DataLazy::collapseToReady()
883      case TRACE:      case TRACE:
884      result=left.trace(m_axis_offset);      result=left.trace(m_axis_offset);
885      break;      break;
886        case SWAP:
887        result=left.swapaxes(m_axis_offset, m_transpose);
888        break;
889        case MINVAL:
890        result=left.minval();
891        break;
892        case MAXVAL:
893        result=left.minval();
894        break;
895      default:      default:
896      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)+".");
897    }    }
# Line 745  DataLazy::collapse() Line 919  DataLazy::collapse()
919    m_op=IDENTITY;    m_op=IDENTITY;
920  }  }
921    
922  /*  
923    \brief Compute the value of the expression (unary operation) for the given sample.  
924    \return Vector which stores the value of the subexpression for the given sample.  
925    \param v A vector to store intermediate results.  
926    \param offset Index in v to begin storing results.  
927    \param sampleNo Sample number to evaluate.  #define PROC_OP(TYPE,X)                               \
928    \param roffset (output parameter) the offset in the return vector where the result begins.      for (int j=0;j<onumsteps;++j)\
929        {\
930    The return value will be an existing vector so do not deallocate it.        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
931    If the result is stored in v it should be stored at the offset given.        { \
932    Everything from offset to the end of v should be considered available for this method to use.  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
933  */  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
934  DataTypes::ValueType*           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
935  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
936             lroffset+=leftstep; \
937             rroffset+=rightstep; \
938          }\
939          lroffset+=oleftstep;\
940          rroffset+=orightstep;\
941        }
942    
943    
944    // The result will be stored in m_samples
945    // The return value is a pointer to the DataVector, offset is the offset within the return value
946    const DataTypes::ValueType*
947    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
948    {
949    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
950        // collapse so we have a 'E' node or an IDENTITY for some other type
951      if (m_readytype!='E' && m_op!=IDENTITY)
952      {
953        collapse();
954      }
955      if (m_op==IDENTITY)  
956      {
957        const ValueType& vec=m_id->getVectorRO();
958        roffset=m_id->getPointOffset(sampleNo, 0);
959    #ifdef LAZY_STACK_PROF
960    int x;
961    if (&x<stackend[omp_get_thread_num()])
962    {
963           stackend[omp_get_thread_num()]=&x;
964    }
965    #endif
966        return &(vec);
967      }
968      if (m_readytype!='E')
969      {
970        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
971      }
972      if (m_sampleids[tid]==sampleNo)
973      {
974        roffset=tid*m_samplesize;
975        return &(m_samples);        // sample is already resolved
976      }
977      m_sampleids[tid]=sampleNo;
978      switch (getOpgroup(m_op))
979      {
980      case G_UNARY:
981      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
982      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
983      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
984      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
985      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
986      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
987      case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
988      default:
989        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
990      }
991    }
992    
993    const DataTypes::ValueType*
994    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
995  {  {
996      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
997      // since we only have one argument we don't need to think about only      // since we only have one argument we don't need to think about only
998      // processing single points.      // processing single points.
999        // we will also know we won't get identity nodes
1000    if (m_readytype!='E')    if (m_readytype!='E')
1001    {    {
1002      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1003    }    }
1004    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);    if (m_op==IDENTITY)
1005    const double* left=&((*vleft)[roffset]);    {
1006    double* result=&(v[offset]);      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1007    roffset=offset;    }
1008      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1009      const double* left=&((*leftres)[roffset]);
1010      roffset=m_samplesize*tid;
1011      double* result=&(m_samples[roffset]);
1012    switch (m_op)    switch (m_op)
1013    {    {
1014      case SIN:        case SIN:  
# Line 880  DataLazy::resolveUnary(ValueType& v, siz Line 1118  DataLazy::resolveUnary(ValueType& v, siz
1118      default:      default:
1119      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1120    }    }
1121    return &v;    return &(m_samples);
1122  }  }
1123    
1124    
1125    const DataTypes::ValueType*
1126    DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)
1127    {
1128        // we assume that any collapsing has been done before we get here
1129        // since we only have one argument we don't need to think about only
1130        // processing single points.
1131        // we will also know we won't get identity nodes
1132      if (m_readytype!='E')
1133      {
1134        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1135      }
1136      if (m_op==IDENTITY)
1137      {
1138        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1139      }
1140      size_t loffset=0;
1141      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1142    
1143      roffset=m_samplesize*tid;
1144      unsigned int ndpps=getNumDPPSample();
1145      unsigned int psize=DataTypes::noValues(getShape());
1146      double* result=&(m_samples[roffset]);
1147      switch (m_op)
1148      {
1149        case MINVAL:
1150        {
1151          for (unsigned int z=0;z<ndpps;++z)
1152          {
1153            FMin op;
1154            *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1155            loffset+=psize;
1156            result++;
1157          }
1158        }
1159        break;
1160        case MAXVAL:
1161        {
1162          for (unsigned int z=0;z<ndpps;++z)
1163          {
1164          FMax op;
1165          *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1166          loffset+=psize;
1167          result++;
1168          }
1169        }
1170        break;
1171        default:
1172        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1173      }
1174      return &(m_samples);
1175    }
1176    
1177    const DataTypes::ValueType*
1178    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
 /*  
   \brief Compute the value of the expression (unary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
1179  {  {
1180      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1181      // since we only have one argument we don't need to think about only      // since we only have one argument we don't need to think about only
1182      // processing single points.      // processing single points.
1183    if (m_readytype!='E')    if (m_readytype!='E')
1184    {    {
1185      throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1186    }    }
1187      // since we can't write the result over the input, we need a result offset further along    if (m_op==IDENTITY)
1188    size_t subroffset=roffset+m_samplesize;    {
1189    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1190    roffset=offset;    }
1191      size_t subroffset;
1192      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1193      roffset=m_samplesize*tid;
1194      size_t loop=0;
1195      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1196      size_t step=getNoValues();
1197      size_t offset=roffset;
1198    switch (m_op)    switch (m_op)
1199    {    {
1200      case SYM:      case SYM:
1201      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1202        {
1203            DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1204            subroffset+=step;
1205            offset+=step;
1206        }
1207      break;      break;
1208      case NSYM:      case NSYM:
1209      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1210        {
1211            DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1212            subroffset+=step;
1213            offset+=step;
1214        }
1215      break;      break;
1216      default:      default:
1217      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1218    }    }
1219    return &v;    return &m_samples;
1220  }  }
1221    
1222  /*  const DataTypes::ValueType*
1223    \brief Compute the value of the expression (unary operation) for the given sample.  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
1224  {  {
1225      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1226      // since we only have one argument we don't need to think about only      // since we only have one argument we don't need to think about only
1227      // processing single points.      // processing single points.
1228    if (m_readytype!='E')    if (m_readytype!='E')
1229    {    {
1230      throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1231    }    }
1232      // since we can't write the result over the input, we need a result offset further along    if (m_op==IDENTITY)
1233    size_t subroffset=roffset+m_samplesize;    {
1234    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);      throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1235    roffset=offset;    }
1236      size_t subroffset;
1237      size_t offset;
1238      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1239      roffset=m_samplesize*tid;
1240      offset=roffset;
1241      size_t loop=0;
1242      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1243      size_t outstep=getNoValues();
1244      size_t instep=m_left->getNoValues();
1245    switch (m_op)    switch (m_op)
1246    {    {
1247      case TRACE:      case TRACE:
1248           DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
1249        {
1250                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1251            subroffset+=instep;
1252            offset+=outstep;
1253        }
1254      break;      break;
1255      case TRANS:      case TRANS:
1256           DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
1257        {
1258                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1259            subroffset+=instep;
1260            offset+=outstep;
1261        }
1262      break;      break;
1263      default:      default:
1264      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1265    }    }
1266    return &v;    return &m_samples;
1267  }  }
1268    
1269    
1270  #define PROC_OP(TYPE,X)                               \  const DataTypes::ValueType*
1271      for (int i=0;i<steps;++i,resultp+=resultStep) \  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1272      { \  {
1273         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    if (m_readytype!='E')
1274         lroffset+=leftStep; \    {
1275         rroffset+=rightStep; \      throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1276      }
1277      if (m_op==IDENTITY)
1278      {
1279        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1280      }
1281      size_t subroffset;
1282      size_t offset;
1283      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1284      roffset=m_samplesize*tid;
1285      offset=roffset;
1286      size_t loop=0;
1287      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1288      size_t outstep=getNoValues();
1289      size_t instep=m_left->getNoValues();
1290      switch (m_op)
1291      {
1292        case SWAP:
1293        for (loop=0;loop<numsteps;++loop)
1294        {
1295                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1296            subroffset+=instep;
1297            offset+=outstep;
1298      }      }
1299        break;
1300        default:
1301        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1302      }
1303      return &m_samples;
1304    }
1305    
1306    
1307    
 /*  
   \brief Compute the value of the expression (binary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
1308  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1309  // have already been collapsed to IDENTITY. So we must have at least one expanded child.  // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1310  // If both children are expanded, then we can process them in a single operation (we treat  // If both children are expanded, then we can process them in a single operation (we treat
# Line 998  DataLazy::resolveNP1OUT_P(ValueType& v, Line 1314  DataLazy::resolveNP1OUT_P(ValueType& v,
1314  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1315  // For example, 2+Vector.  // For example, 2+Vector.
1316  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1317  DataTypes::ValueType*  const DataTypes::ValueType*
1318  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1319  {  {
1320  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1321    
# Line 1013  LAZYDEBUG(cout << "Resolve binary: " << Line 1329  LAZYDEBUG(cout << "Resolve binary: " <<
1329    }    }
1330    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1331    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1332    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1333    int steps=(bigloops?1:getNumDPPSample());    {
1334    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint      throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1335    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1336    {    size_t leftsize=m_left->getNoValues();
1337      if (!leftScalar && !rightScalar)    size_t rightsize=m_right->getNoValues();
1338      {    size_t chunksize=1;           // how many doubles will be processed in one go
1339         throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");    int leftstep=0;       // how far should the left offset advance after each step
1340      }    int rightstep=0;
1341      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    int numsteps=0;       // total number of steps for the inner loop
1342      chunksize=1;    // for scalar    int oleftstep=0;  // the o variables refer to the outer loop
1343    }        int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1344    int leftStep=((leftExp && (!rightExp || rightScalar))? m_right->getNoValues() : 0);    int onumsteps=1;
1345    int rightStep=((rightExp && (!leftExp || leftScalar))? m_left->getNoValues() : 0);    
1346    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1347      bool RES=(rightExp && rightScalar);
1348      bool LS=(!leftExp && leftScalar); // left is a single scalar
1349      bool RS=(!rightExp && rightScalar);
1350      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1351      bool RN=(!rightExp && !rightScalar);
1352      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1353      bool REN=(rightExp && !rightScalar);
1354    
1355      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1356      {
1357        chunksize=m_left->getNumDPPSample()*leftsize;
1358        leftstep=0;
1359        rightstep=0;
1360        numsteps=1;
1361      }
1362      else if (LES || RES)
1363      {
1364        chunksize=1;
1365        if (LES)        // left is an expanded scalar
1366        {
1367            if (RS)
1368            {
1369               leftstep=1;
1370               rightstep=0;
1371               numsteps=m_left->getNumDPPSample();
1372            }
1373            else        // RN or REN
1374            {
1375               leftstep=0;
1376               oleftstep=1;
1377               rightstep=1;
1378               orightstep=(RN ? -(int)rightsize : 0);
1379               numsteps=rightsize;
1380               onumsteps=m_left->getNumDPPSample();
1381            }
1382        }
1383        else        // right is an expanded scalar
1384        {
1385            if (LS)
1386            {
1387               rightstep=1;
1388               leftstep=0;
1389               numsteps=m_right->getNumDPPSample();
1390            }
1391            else
1392            {
1393               rightstep=0;
1394               orightstep=1;
1395               leftstep=1;
1396               oleftstep=(LN ? -(int)leftsize : 0);
1397               numsteps=leftsize;
1398               onumsteps=m_right->getNumDPPSample();
1399            }
1400        }
1401      }
1402      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1403      {
1404        if (LEN)    // and Right will be a single value
1405        {
1406            chunksize=rightsize;
1407            leftstep=rightsize;
1408            rightstep=0;
1409            numsteps=m_left->getNumDPPSample();
1410            if (RS)
1411            {
1412               numsteps*=leftsize;
1413            }
1414        }
1415        else    // REN
1416        {
1417            chunksize=leftsize;
1418            rightstep=leftsize;
1419            leftstep=0;
1420            numsteps=m_right->getNumDPPSample();
1421            if (LS)
1422            {
1423               numsteps*=rightsize;
1424            }
1425        }
1426      }
1427    
1428      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1429      // Get the values of sub-expressions      // Get the values of sub-expressions
1430    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
1431    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1432      // the right child starts further along.  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1433    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1434    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1435    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1436    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1437    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1438    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1439    
1440    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1441    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1442    
1443    
1444      roffset=m_samplesize*tid;
1445      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
1446    switch(m_op)    switch(m_op)
1447    {    {
1448      case ADD:      case ADD:
# Line 1053  LAZYDEBUG(cout << "Resolve binary: " << Line 1463  LAZYDEBUG(cout << "Resolve binary: " <<
1463      default:      default:
1464      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1465    }    }
1466    roffset=offset;    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1467    return &v;    return &m_samples;
1468  }  }
1469    
1470    
 /*  
   \brief Compute the value of the expression (tensor product) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
1471  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1472  // have already been collapsed to IDENTITY. So we must have at least one expanded child.  // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1473  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1474  DataTypes::ValueType*  const DataTypes::ValueType*
1475  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
1476  {  {
1477  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1478    
# Line 1083  LAZYDEBUG(cout << "Resolve TensorProduct Line 1481  LAZYDEBUG(cout << "Resolve TensorProduct
1481    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1482    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1483    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1484    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1485    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1486    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0  
1487      // Get the values of sub-expressions (leave a gap of one sample for the result).    int resultStep=getNoValues();
1488    const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);    roffset=m_samplesize*tid;
1489    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);    size_t offset=roffset;
1490    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved  
1491      const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1492    
1493      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1494    
1495    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1496    cout << getNoValues() << endl;)
1497    
1498    
1499    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1500    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1501    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1502    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1503    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1504    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1505    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1506    
1507      double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
1508    switch(m_op)    switch(m_op)
1509    {    {
1510      case PROD:      case PROD:
# Line 1097  LAZYDEBUG(cout << "Resolve TensorProduct Line 1512  LAZYDEBUG(cout << "Resolve TensorProduct
1512      {      {
1513            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1514            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1515    
1516    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1517    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1518    
1519            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);
1520    
1521        lroffset+=leftStep;        lroffset+=leftStep;
1522        rroffset+=rightStep;        rroffset+=rightStep;
1523      }      }
# Line 1106  LAZYDEBUG(cout << "Resolve TensorProduct Line 1526  LAZYDEBUG(cout << "Resolve TensorProduct
1526      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1527    }    }
1528    roffset=offset;    roffset=offset;
1529    return &v;    return &m_samples;
1530  }  }
1531    
1532    
1533    const DataTypes::ValueType*
1534    DataLazy::resolveSample(int sampleNo, size_t& roffset)
1535    {
1536    #ifdef _OPENMP
1537        int tid=omp_get_thread_num();
1538    #else
1539        int tid=0;
1540    #endif
1541    
1542  /*  #ifdef LAZY_STACK_PROF
1543    \brief Compute the value of the expression for the given sample.      stackstart[tid]=&tid;
1544    \return Vector which stores the value of the subexpression for the given sample.      stackend[tid]=&tid;
1545    \param v A vector to store intermediate results.      const DataTypes::ValueType* r=resolveNodeSample(tid, sampleNo, roffset);
1546    \param offset Index in v to begin storing results.      size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1547    \param sampleNo Sample number to evaluate.      #pragma omp critical
1548    \param roffset (output parameter) the offset in the return vector where the result begins.      if (d>maxstackuse)
1549        {
1550    cout << "Max resolve Stack use " << d << endl;
1551            maxstackuse=d;
1552        }
1553        return r;
1554    #else
1555        return resolveNodeSample(tid, sampleNo, roffset);
1556    #endif
1557    }
1558    
   The return value will be an existing vector so do not deallocate it.  
 */  
 // the vector and the offset are a place where the method could write its data if it wishes  
 // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
 // Hence the return value to indicate where the data is actually stored.  
 // Regardless, the storage should be assumed to be used, even if it isn't.  
1559    
1560  // the roffset is the offset within the returned vector where the data begins  // This needs to do the work of the identity constructor
1561  const DataTypes::ValueType*  void
1562  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveToIdentity()
1563  {  {
1564  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)     if (m_op==IDENTITY)
1565      // collapse so we have a 'E' node or an IDENTITY for some other type      return;
1566    if (m_readytype!='E' && m_op!=IDENTITY)     DataReady_ptr p=resolveNodeWorker();
1567    {     makeIdentity(p);
1568      collapse();  }
1569    }  
1570    if (m_op==IDENTITY)    void DataLazy::makeIdentity(const DataReady_ptr& p)
1571    {  {
1572      const ValueType& vec=m_id->getVector();     m_op=IDENTITY;
1573      if (m_readytype=='C')     m_axis_offset=0;
1574      {     m_transpose=0;
1575      roffset=0;     m_SL=m_SM=m_SR=0;
1576      return &(vec);     m_children=m_height=0;
1577      }     m_id=p;
1578      roffset=m_id->getPointOffset(sampleNo, 0);     if(p->isConstant()) {m_readytype='C';}
1579      return &(vec);     else if(p->isExpanded()) {m_readytype='E';}
1580    }     else if (p->isTagged()) {m_readytype='T';}
1581    if (m_readytype!='E')     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1582    {     m_samplesize=p->getNumDPPSample()*p->getNoValues();
1583      throw DataException("Programmer Error - Collapse did not produce an expanded node.");     m_left.reset();
1584    }     m_right.reset();
   switch (getOpgroup(m_op))  
   {  
   case G_UNARY:  
   case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);  
   case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);  
   case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);  
   case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);  
   case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);  
   default:  
     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
   }  
1585  }  }
1586    
1587    
 // To simplify the memory management, all threads operate on one large vector, rather than one each.  
 // Each sample is evaluated independently and copied into the result DataExpanded.  
1588  DataReady_ptr  DataReady_ptr
1589  DataLazy::resolve()  DataLazy::resolve()
1590  {  {
1591        resolveToIdentity();
1592        return m_id;
1593    }
1594    
1595  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)  // This version of resolve uses storage in each node to hold results
1596  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)  DataReady_ptr
1597    DataLazy::resolveNodeWorker()
1598    {
1599    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
1600    {    {
1601      collapse();      collapse();
# Line 1183  LAZYDEBUG(cout << "Buffers=" << m_buffsR Line 1605  LAZYDEBUG(cout << "Buffers=" << m_buffsR
1605      return m_id;      return m_id;
1606    }    }
1607      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
   size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough  
     // storage to evaluate its expression  
   int numthreads=1;  
 #ifdef _OPENMP  
   numthreads=getNumberOfThreads();  
 #endif  
   ValueType v(numthreads*threadbuffersize);  
 LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)  
1608    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1609    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
1610    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1611    
1612    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1613    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1614    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Storage for answer
1615    size_t resoffset=0;       // where in the vector to find the answer  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1616    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel private(sample,res)
1617    for (sample=0;sample<totalsamples;++sample)    {
1618    {      size_t roffset=0;
1619  LAZYDEBUG(cout << "################################# " << sample << endl;)  #ifdef LAZY_STACK_PROF
1620        stackstart[omp_get_thread_num()]=&roffset;
1621        stackend[omp_get_thread_num()]=&roffset;
1622    #endif
1623        #pragma omp for schedule(static)
1624        for (sample=0;sample<totalsamples;++sample)
1625        {
1626            roffset=0;
1627  #ifdef _OPENMP  #ifdef _OPENMP
1628      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);              res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1629  #else  #else
1630      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.              res=resolveNodeSample(0,sample,roffset);
1631  #endif  #endif
1632  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1633      outoffset=result->getPointOffset(sample,0);  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1634  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)              DataVector::size_type outoffset=result->getPointOffset(sample,0);
1635      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector              memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
1636      {      }
1637      resvec[outoffset]=(*res)[resoffset];    }
1638      }  #ifdef LAZY_STACK_PROF
1639  LAZYDEBUG(cerr << "*********************************" << endl;)    for (int i=0;i<getNumberOfThreads();++i)
1640      {
1641        size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1642    //  cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1643        if (r>maxstackuse)
1644        {
1645            maxstackuse=r;
1646        }
1647    }    }
1648      cout << "Max resolve Stack use=" << maxstackuse << endl;
1649    #endif
1650    return resptr;    return resptr;
1651  }  }
1652    
# Line 1224  std::string Line 1654  std::string
1654  DataLazy::toString() const  DataLazy::toString() const
1655  {  {
1656    ostringstream oss;    ostringstream oss;
1657    oss << "Lazy Data:";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1658    intoString(oss);    if (escriptParams.getPRINT_LAZY_TREE()==0)
1659      {
1660          intoString(oss);
1661      }
1662      else
1663      {
1664        oss << endl;
1665        intoTreeString(oss,"");
1666      }
1667    return oss.str();    return oss.str();
1668  }  }
1669    
# Line 1233  DataLazy::toString() const Line 1671  DataLazy::toString() const
1671  void  void
1672  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1673  {  {
1674    //    oss << "[" << m_children <<";"<<m_height <<"]";
1675    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1676    {    {
1677    case G_IDENTITY:    case G_IDENTITY:
# Line 1265  DataLazy::intoString(ostringstream& oss) Line 1704  DataLazy::intoString(ostringstream& oss)
1704    case G_UNARY_P:    case G_UNARY_P:
1705    case G_NP1OUT:    case G_NP1OUT:
1706    case G_NP1OUT_P:    case G_NP1OUT_P:
1707      case G_REDUCTION:
1708      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1709      m_left->intoString(oss);      m_left->intoString(oss);
1710      oss << ')';      oss << ')';
# Line 1276  DataLazy::intoString(ostringstream& oss) Line 1716  DataLazy::intoString(ostringstream& oss)
1716      m_right->intoString(oss);      m_right->intoString(oss);
1717      oss << ')';      oss << ')';
1718      break;      break;
1719      case G_NP1OUT_2P:
1720        oss << opToString(m_op) << '(';
1721        m_left->intoString(oss);
1722        oss << ", " << m_axis_offset << ", " << m_transpose;
1723        oss << ')';
1724        break;
1725    default:    default:
1726      oss << "UNKNOWN";      oss << "UNKNOWN";
1727    }    }
1728  }  }
1729    
1730    
1731    void
1732    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1733    {
1734      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1735      switch (getOpgroup(m_op))
1736      {
1737      case G_IDENTITY:
1738        if (m_id->isExpanded())
1739        {
1740           oss << "E";
1741        }
1742        else if (m_id->isTagged())
1743        {
1744          oss << "T";
1745        }
1746        else if (m_id->isConstant())
1747        {
1748          oss << "C";
1749        }
1750        else
1751        {
1752          oss << "?";
1753        }
1754        oss << '@' << m_id.get() << endl;
1755        break;
1756      case G_BINARY:
1757        oss << opToString(m_op) << endl;
1758        indent+='.';
1759        m_left->intoTreeString(oss, indent);
1760        m_right->intoTreeString(oss, indent);
1761        break;
1762      case G_UNARY:
1763      case G_UNARY_P:
1764      case G_NP1OUT:
1765      case G_NP1OUT_P:
1766      case G_REDUCTION:
1767        oss << opToString(m_op) << endl;
1768        indent+='.';
1769        m_left->intoTreeString(oss, indent);
1770        break;
1771      case G_TENSORPROD:
1772        oss << opToString(m_op) << endl;
1773        indent+='.';
1774        m_left->intoTreeString(oss, indent);
1775        m_right->intoTreeString(oss, indent);
1776        break;
1777      case G_NP1OUT_2P:
1778        oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1779        indent+='.';
1780        m_left->intoTreeString(oss, indent);
1781        break;
1782      default:
1783        oss << "UNKNOWN";
1784      }
1785    }
1786    
1787    
1788  DataAbstract*  DataAbstract*
1789  DataLazy::deepCopy()  DataLazy::deepCopy()
1790  {  {
1791    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1792    {    {
1793    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1794    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY:
1795      case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1796      case G_UNARY_P:   return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1797    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);    case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1798    case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);    case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1799    case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);    case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1800      case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
1801      case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1802    default:    default:
1803      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");      throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1804    }    }
1805  }  }
1806    
1807    
1808    
1809  // There is no single, natural interpretation of getLength on DataLazy.  // There is no single, natural interpretation of getLength on DataLazy.
1810  // Instances of DataReady can look at the size of their vectors.  // Instances of DataReady can look at the size of their vectors.
1811  // For lazy though, it could be the size the data would be if it were resolved;  // For lazy though, it could be the size the data would be if it were resolved;
# Line 1372  DataLazy::getPointOffset(int sampleNo, Line 1881  DataLazy::getPointOffset(int sampleNo,
1881    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).");
1882  }  }
1883    
1884  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
1885  // 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.  
1886  void  void
1887  DataLazy::setToZero()  DataLazy::setToZero()
1888  {  {
1889    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
1890    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1891    m_op=IDENTITY;  //   m_op=IDENTITY;
1892    m_right.reset();    //   m_right.reset();  
1893    m_left.reset();  //   m_left.reset();
1894    m_readytype='C';  //   m_readytype='C';
1895    m_buffsRequired=1;  //   m_buffsRequired=1;
1896    
1897      privdebug=privdebug;  // to stop the compiler complaining about unused privdebug
1898      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1899    }
1900    
1901    bool
1902    DataLazy::actsExpanded() const
1903    {
1904        return (m_readytype=='E');
1905  }  }
1906    
1907  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26