/[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 2770 by jfenwick, Wed Nov 25 01:24:51 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    #include "EscriptParams.h"
32    
33    #include <iomanip>      // for some fancy formatting in debug
34    
35    // #define LAZYDEBUG(X) if (privdebug){X;}
36    #define LAZYDEBUG(X)
37    namespace
38    {
39    bool privdebug=false;
40    
41    #define ENABLEDEBUG privdebug=true;
42    #define DISABLEDEBUG privdebug=false;
43    }
44    
45    // #define SIZELIMIT 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?
52  ~~~~~~~~~~~~~~~~~~~~~~~  ~~~~~~~~~~~~~~~~~~~~~~~
# Line 70  The convention that I use, is that the r Line 89  The convention that I use, is that the r
89  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.
90  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.
91    
92  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):
93  1) Add to the ES_optype.  1) Add to the ES_optype.
94  2) determine what opgroup your operation belongs to (X)  2) determine what opgroup your operation belongs to (X)
95  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 115  enum ES_opgroup
115     G_IDENTITY,     G_IDENTITY,
116     G_BINARY,        // pointwise operations with two arguments     G_BINARY,        // pointwise operations with two arguments
117     G_UNARY,     // pointwise operations with one argument     G_UNARY,     // pointwise operations with one argument
118       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
119     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
120     G_TENSORPROD     // general tensor product     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
121       G_TENSORPROD,    // general tensor product
122       G_NP1OUT_2P,     // non-pointwise op with one output requiring two params
123       G_REDUCTION      // non-pointwise unary op with a scalar output
124  };  };
125    
126    
# Line 108  string ES_opstrings[]={"UNKNOWN","IDENTI Line 131  string ES_opstrings[]={"UNKNOWN","IDENTI
131              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
132              "asinh","acosh","atanh",              "asinh","acosh","atanh",
133              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
134              "1/","where>0","where<0","where>=0","where<=0",              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
135              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
136              "prod"};              "prod",
137  int ES_opcount=36;              "transpose", "trace",
138                "swapaxes",
139                "minval", "maxval"};
140    int ES_opcount=43;
141  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,
142              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
143              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
144              G_UNARY,G_UNARY,G_UNARY,                    // 20              G_UNARY,G_UNARY,G_UNARY,                    // 20
145              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
146              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
147              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
148              G_TENSORPROD};              G_TENSORPROD,
149                G_NP1OUT_P, G_NP1OUT_P,
150                G_NP1OUT_2P,
151                G_REDUCTION, G_REDUCTION};
152  inline  inline
153  ES_opgroup  ES_opgroup
154  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 164  resultShape(DataAbstract_ptr left, DataA Line 193  resultShape(DataAbstract_ptr left, DataA
193        {        {
194          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.");
195        }        }
196    
197        if (left->getRank()==0)   // we need to allow scalar * anything        if (left->getRank()==0)   // we need to allow scalar * anything
198        {        {
199          return right->getShape();          return right->getShape();
# Line 177  resultShape(DataAbstract_ptr left, DataA Line 207  resultShape(DataAbstract_ptr left, DataA
207      return left->getShape();      return left->getShape();
208  }  }
209    
210    // return the shape for "op left"
211    
212    DataTypes::ShapeType
213    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
214    {
215        switch(op)
216        {
217            case TRANS:
218           {            // for the scoping of variables
219            const DataTypes::ShapeType& s=left->getShape();
220            DataTypes::ShapeType sh;
221            int rank=left->getRank();
222            if (axis_offset<0 || axis_offset>rank)
223            {
224                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
225                }
226                for (int i=0; i<rank; i++)
227            {
228               int index = (axis_offset+i)%rank;
229                   sh.push_back(s[index]); // Append to new shape
230                }
231            return sh;
232           }
233        break;
234        case TRACE:
235           {
236            int rank=left->getRank();
237            if (rank<2)
238            {
239               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
240            }
241            if ((axis_offset>rank-2) || (axis_offset<0))
242            {
243               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
244            }
245            if (rank==2)
246            {
247               return DataTypes::scalarShape;
248            }
249            else if (rank==3)
250            {
251               DataTypes::ShapeType sh;
252                   if (axis_offset==0)
253               {
254                    sh.push_back(left->getShape()[2]);
255                   }
256                   else     // offset==1
257               {
258                sh.push_back(left->getShape()[0]);
259                   }
260               return sh;
261            }
262            else if (rank==4)
263            {
264               DataTypes::ShapeType sh;
265               const DataTypes::ShapeType& s=left->getShape();
266                   if (axis_offset==0)
267               {
268                    sh.push_back(s[2]);
269                    sh.push_back(s[3]);
270                   }
271                   else if (axis_offset==1)
272               {
273                    sh.push_back(s[0]);
274                    sh.push_back(s[3]);
275                   }
276               else     // offset==2
277               {
278                sh.push_back(s[0]);
279                sh.push_back(s[1]);
280               }
281               return sh;
282            }
283            else        // unknown rank
284            {
285               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
286            }
287           }
288        break;
289            default:
290        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
291        }
292    }
293    
294    DataTypes::ShapeType
295    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
296    {
297         // This code taken from the Data.cpp swapaxes() method
298         // Some of the checks are probably redundant here
299         int axis0_tmp,axis1_tmp;
300         const DataTypes::ShapeType& s=left->getShape();
301         DataTypes::ShapeType out_shape;
302         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
303         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
304         int rank=left->getRank();
305         if (rank<2) {
306            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
307         }
308         if (axis0<0 || axis0>rank-1) {
309            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
310         }
311         if (axis1<0 || axis1>rank-1) {
312             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
313         }
314         if (axis0 == axis1) {
315             throw DataException("Error - Data::swapaxes: axis indices must be different.");
316         }
317         if (axis0 > axis1) {
318             axis0_tmp=axis1;
319             axis1_tmp=axis0;
320         } else {
321             axis0_tmp=axis0;
322             axis1_tmp=axis1;
323         }
324         for (int i=0; i<rank; i++) {
325           if (i == axis0_tmp) {
326              out_shape.push_back(s[axis1_tmp]);
327           } else if (i == axis1_tmp) {
328              out_shape.push_back(s[axis0_tmp]);
329           } else {
330              out_shape.push_back(s[i]);
331           }
332         }
333        return out_shape;
334    }
335    
336    
337  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
338  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
339  // 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 354  GTPShape(DataAbstract_ptr left, DataAbst
354    else if (transpose == 2)  { start1 = rank1-axis_offset; }    else if (transpose == 2)  { start1 = rank1-axis_offset; }
355    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"); }
356    
357      if (rank0<axis_offset)
358      {
359        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
360      }
361    
362    // Adjust the shapes for transpose    // Adjust the shapes for transpose
363    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 387  GTPShape(DataAbstract_ptr left, DataAbst
387       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
388       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
389    }    }
   return shape2;  
 }  
390    
391      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
392      {
393         ostringstream os;
394         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
395         throw DataException(os.str());
396      }
397    
398  // determine the number of points in the result of "left op right"    return shape2;
 // 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: return max(left->getBuffsRequired(),1);  
    case G_NP1OUT: 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)+".");  
    }  
399  }  }
400    
   
401  }   // end anonymous namespace  }   // end anonymous namespace
402    
403    
# Line 279  opToString(ES_optype op) Line 413  opToString(ES_optype op)
413    return ES_opstrings[op];    return ES_opstrings[op];
414  }  }
415    
416    void DataLazy::LazyNodeSetup()
417    {
418    #ifdef _OPENMP
419        int numthreads=omp_get_max_threads();
420        m_samples.resize(numthreads*m_samplesize);
421        m_sampleids=new int[numthreads];
422        for (int i=0;i<numthreads;++i)
423        {
424            m_sampleids[i]=-1;  
425        }
426    #else
427        m_samples.resize(m_samplesize);
428        m_sampleids=new int[1];
429        m_sampleids[0]=-1;
430    #endif  // _OPENMP
431    }
432    
433    
434    // Creates an identity node
435  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
436      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
437      m_op(IDENTITY),      ,m_sampleids(0),
438      m_axis_offset(0),      m_samples(1)
     m_transpose(0),  
     m_SL(0), m_SM(0), m_SR(0)  
439  {  {
440     if (p->isLazy())     if (p->isLazy())
441     {     {
# Line 296  DataLazy::DataLazy(DataAbstract_ptr p) Line 446  DataLazy::DataLazy(DataAbstract_ptr p)
446     }     }
447     else     else
448     {     {
449      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
450      if(p->isConstant()) {m_readytype='C';}      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
451      else if(p->isExpanded()) {m_readytype='E';}      makeIdentity(dr);
452      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.");}  
453     }     }
454     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;  
455  }  }
456    
   
   
   
457  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
458      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
459      m_op(op),      m_op(op),
460      m_axis_offset(0),      m_axis_offset(0),
461      m_transpose(0),      m_transpose(0),
462      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
463  {  {
464     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
465     {     {
466      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.");
467     }     }
# Line 334  DataLazy::DataLazy(DataAbstract_ptr left Line 477  DataLazy::DataLazy(DataAbstract_ptr left
477     }     }
478     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
479     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
480     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
481     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
482       m_height=m_left->m_height+1;
483       LazyNodeSetup();
484       SIZELIMIT
485  }  }
486    
487    
# Line 346  DataLazy::DataLazy(DataAbstract_ptr left Line 491  DataLazy::DataLazy(DataAbstract_ptr left
491      m_op(op),      m_op(op),
492      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
493  {  {
494    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
495     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
496     {     {
497      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 508  DataLazy::DataLazy(DataAbstract_ptr left
508     {     {
509      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
510      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
511    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
512     }     }
513     left->operandCheck(*right);     left->operandCheck(*right);
514    
515     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
516     {     {
517      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
518    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
519     }     }
520     else     else
521     {     {
522      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
523    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
524     }     }
525     if (right->isLazy())     if (right->isLazy())
526     {     {
527      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
528    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
529     }     }
530     else     else
531     {     {
532      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
533    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
534     }     }
535     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
536     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 396  DataLazy::DataLazy(DataAbstract_ptr left Line 547  DataLazy::DataLazy(DataAbstract_ptr left
547      m_readytype='C';      m_readytype='C';
548     }     }
549     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
550     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_children=m_left->m_children+m_right->m_children+2;
551     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
552  cout << "(3)Lazy created with " << m_samplesize << endl;     LazyNodeSetup();
553       SIZELIMIT
554    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
555  }  }
556    
557  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 580  DataLazy::DataLazy(DataAbstract_ptr left
580      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
581      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
582     }     }
583     left->operandCheck(*right);  //    left->operandCheck(*right);
584    
585     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
586     {     {
# Line 460  DataLazy::DataLazy(DataAbstract_ptr left Line 613  DataLazy::DataLazy(DataAbstract_ptr left
613      m_readytype='C';      m_readytype='C';
614     }     }
615     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
616     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_children=m_left->m_children+m_right->m_children+2;
617     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
618  cout << "(4)Lazy created with " << m_samplesize << endl;     LazyNodeSetup();
619       SIZELIMIT
620    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
621  }  }
622    
623    
624  DataLazy::~DataLazy()  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
625        : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
626        m_op(op),
627        m_axis_offset(axis_offset),
628        m_transpose(0),
629        m_tol(0)
630  {  {
631       if ((getOpgroup(op)!=G_NP1OUT_P))
632       {
633        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
634       }
635       DataLazy_ptr lleft;
636       if (!left->isLazy())
637       {
638        lleft=DataLazy_ptr(new DataLazy(left));
639       }
640       else
641       {
642        lleft=dynamic_pointer_cast<DataLazy>(left);
643       }
644       m_readytype=lleft->m_readytype;
645       m_left=lleft;
646       m_samplesize=getNumDPPSample()*getNoValues();
647       m_children=m_left->m_children+1;
648       m_height=m_left->m_height+1;
649       LazyNodeSetup();
650       SIZELIMIT
651    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
652  }  }
653    
654    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
655  int      : parent(left->getFunctionSpace(), left->getShape()),
656  DataLazy::getBuffsRequired() const      m_op(op),
657        m_axis_offset(0),
658        m_transpose(0),
659        m_tol(tol)
660  {  {
661      return m_buffsRequired;     if ((getOpgroup(op)!=G_UNARY_P))
662       {
663        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
664       }
665       DataLazy_ptr lleft;
666       if (!left->isLazy())
667       {
668        lleft=DataLazy_ptr(new DataLazy(left));
669       }
670       else
671       {
672        lleft=dynamic_pointer_cast<DataLazy>(left);
673       }
674       m_readytype=lleft->m_readytype;
675       m_left=lleft;
676       m_samplesize=getNumDPPSample()*getNoValues();
677       m_children=m_left->m_children+1;
678       m_height=m_left->m_height+1;
679       LazyNodeSetup();
680       SIZELIMIT
681    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
682  }  }
683    
684    
685  size_t  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
686  DataLazy::getMaxSampleSize() const      : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
687        m_op(op),
688        m_axis_offset(axis0),
689        m_transpose(axis1),
690        m_tol(0)
691    {
692       if ((getOpgroup(op)!=G_NP1OUT_2P))
693       {
694        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
695       }
696       DataLazy_ptr lleft;
697       if (!left->isLazy())
698       {
699        lleft=DataLazy_ptr(new DataLazy(left));
700       }
701       else
702       {
703        lleft=dynamic_pointer_cast<DataLazy>(left);
704       }
705       m_readytype=lleft->m_readytype;
706       m_left=lleft;
707       m_samplesize=getNumDPPSample()*getNoValues();
708       m_children=m_left->m_children+1;
709       m_height=m_left->m_height+1;
710       LazyNodeSetup();
711       SIZELIMIT
712    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
713    }
714    
715    DataLazy::~DataLazy()
716  {  {
717      return m_maxsamplesize;     delete[] m_sampleids;
718  }  }
719    
720    
721  /*  /*
722    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
723    This does the work for the collapse method.    This does the work for the collapse method.
# Line 602  DataLazy::collapseToReady() Line 836  DataLazy::collapseToReady()
836      case LEZ:      case LEZ:
837      result=left.whereNonPositive();      result=left.whereNonPositive();
838      break;      break;
839        case NEZ:
840        result=left.whereNonZero(m_tol);
841        break;
842        case EZ:
843        result=left.whereZero(m_tol);
844        break;
845      case SYM:      case SYM:
846      result=left.symmetric();      result=left.symmetric();
847      break;      break;
# Line 611  DataLazy::collapseToReady() Line 851  DataLazy::collapseToReady()
851      case PROD:      case PROD:
852      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);      result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
853      break;      break;
854        case TRANS:
855        result=left.transpose(m_axis_offset);
856        break;
857        case TRACE:
858        result=left.trace(m_axis_offset);
859        break;
860        case SWAP:
861        result=left.swapaxes(m_axis_offset, m_transpose);
862        break;
863        case MINVAL:
864        result=left.minval();
865        break;
866        case MAXVAL:
867        result=left.minval();
868        break;
869      default:      default:
870      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)+".");
871    }    }
# Line 638  DataLazy::collapse() Line 893  DataLazy::collapse()
893    m_op=IDENTITY;    m_op=IDENTITY;
894  }  }
895    
896  /*  
897    \brief Compute the value of the expression (unary operation) for the given sample.  
898    \return Vector which stores the value of the subexpression for the given sample.  
899    \param v A vector to store intermediate results.  
900    \param offset Index in v to begin storing results.  
901    \param sampleNo Sample number to evaluate.  #define PROC_OP(TYPE,X)                               \
902    \param roffset (output parameter) the offset in the return vector where the result begins.      for (int j=0;j<onumsteps;++j)\
903        {\
904    The return value will be an existing vector so do not deallocate it.        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
905    If the result is stored in v it should be stored at the offset given.        { \
906    Everything from offset to the end of v should be considered available for this method to use.  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
907  */  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
908  DataTypes::ValueType*           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
909  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
910             lroffset+=leftstep; \
911             rroffset+=rightstep; \
912          }\
913          lroffset+=oleftstep;\
914          rroffset+=orightstep;\
915        }
916    
917    
918    // The result will be stored in m_samples
919    // The return value is a pointer to the DataVector, offset is the offset within the return value
920    const DataTypes::ValueType*
921    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
922    {
923    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
924        // collapse so we have a 'E' node or an IDENTITY for some other type
925      if (m_readytype!='E' && m_op!=IDENTITY)
926      {
927        collapse();
928      }
929      if (m_op==IDENTITY)  
930      {
931        const ValueType& vec=m_id->getVectorRO();
932        roffset=m_id->getPointOffset(sampleNo, 0);
933        return &(vec);
934      }
935      if (m_readytype!='E')
936      {
937        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
938      }
939      if (m_sampleids[tid]==sampleNo)
940      {
941        roffset=tid*m_samplesize;
942        return &(m_samples);        // sample is already resolved
943      }
944      m_sampleids[tid]=sampleNo;
945      switch (getOpgroup(m_op))
946      {
947      case G_UNARY:
948      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
949      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
950      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
951      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
952      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
953      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
954      case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
955      default:
956        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
957      }
958    }
959    
960    const DataTypes::ValueType*
961    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
962  {  {
963      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
964      // 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
965      // processing single points.      // processing single points.
966        // we will also know we won't get identity nodes
967    if (m_readytype!='E')    if (m_readytype!='E')
968    {    {
969      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
970    }    }
971    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);    if (m_op==IDENTITY)
972    const double* left=&((*vleft)[roffset]);    {
973    double* result=&(v[offset]);      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
974    roffset=offset;    }
975      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
976      const double* left=&((*leftres)[roffset]);
977      roffset=m_samplesize*tid;
978      double* result=&(m_samples[roffset]);
979    switch (m_op)    switch (m_op)
980    {    {
981      case SIN:        case SIN:  
# Line 762  DataLazy::resolveUnary(ValueType& v, siz Line 1074  DataLazy::resolveUnary(ValueType& v, siz
1074      case LEZ:      case LEZ:
1075      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));
1076      break;      break;
1077    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1078        case NEZ:
1079        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1080        break;
1081        case EZ:
1082        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1083        break;
1084    
1085      default:      default:
1086      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1087    }    }
1088    return &v;    return &(m_samples);
1089  }  }
1090    
1091    
1092  /*  const DataTypes::ValueType*
1093    \brief Compute the value of the expression (unary operation) for the given sample.  DataLazy::resolveNodeReduction(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(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
1094  {  {
1095      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1096      // 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
1097      // processing single points.      // processing single points.
1098        // we will also know we won't get identity nodes
1099    if (m_readytype!='E')    if (m_readytype!='E')
1100    {    {
1101      throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1102    }    }
1103      // since we can't write the result over the input, we need a result offset further along    if (m_op==IDENTITY)
1104    size_t subroffset=roffset+m_samplesize;    {
1105    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1106    roffset=offset;    }
1107      size_t loffset=0;
1108      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1109    
1110      roffset=m_samplesize*tid;
1111      unsigned int ndpps=getNumDPPSample();
1112      unsigned int psize=DataTypes::noValues(getShape());
1113      double* result=&(m_samples[roffset]);
1114      switch (m_op)
1115      {
1116        case MINVAL:
1117        {
1118          for (unsigned int z=0;z<ndpps;++z)
1119          {
1120            FMin op;
1121            *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1122            loffset+=psize;
1123            result++;
1124          }
1125        }
1126        break;
1127        case MAXVAL:
1128        {
1129          for (unsigned int z=0;z<ndpps;++z)
1130          {
1131          FMax op;
1132          *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1133          loffset+=psize;
1134          result++;
1135          }
1136        }
1137        break;
1138        default:
1139        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1140      }
1141      return &(m_samples);
1142    }
1143    
1144    const DataTypes::ValueType*
1145    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
1146    {
1147        // we assume that any collapsing has been done before we get here
1148        // since we only have one argument we don't need to think about only
1149        // processing single points.
1150      if (m_readytype!='E')
1151      {
1152        throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1153      }
1154      if (m_op==IDENTITY)
1155      {
1156        throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1157      }
1158      size_t subroffset;
1159      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1160      roffset=m_samplesize*tid;
1161      size_t loop=0;
1162      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1163      size_t step=getNoValues();
1164      size_t offset=roffset;
1165    switch (m_op)    switch (m_op)
1166    {    {
1167      case SYM:      case SYM:
1168      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1169        {
1170            DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1171            subroffset+=step;
1172            offset+=step;
1173        }
1174      break;      break;
1175      case NSYM:      case NSYM:
1176      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1177        {
1178            DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1179            subroffset+=step;
1180            offset+=step;
1181        }
1182      break;      break;
1183      default:      default:
1184      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1185    }    }
1186    return &v;    return &m_samples;
1187  }  }
1188    
1189    const DataTypes::ValueType*
1190    DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
1191    {
1192        // we assume that any collapsing has been done before we get here
1193        // since we only have one argument we don't need to think about only
1194        // processing single points.
1195      if (m_readytype!='E')
1196      {
1197        throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1198      }
1199      if (m_op==IDENTITY)
1200      {
1201        throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1202      }
1203      size_t subroffset;
1204      size_t offset;
1205      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1206      roffset=m_samplesize*tid;
1207      offset=roffset;
1208      size_t loop=0;
1209      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1210      size_t outstep=getNoValues();
1211      size_t instep=m_left->getNoValues();
1212      switch (m_op)
1213      {
1214        case TRACE:
1215        for (loop=0;loop<numsteps;++loop)
1216        {
1217                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1218            subroffset+=instep;
1219            offset+=outstep;
1220        }
1221        break;
1222        case TRANS:
1223        for (loop=0;loop<numsteps;++loop)
1224        {
1225                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1226            subroffset+=instep;
1227            offset+=outstep;
1228        }
1229        break;
1230        default:
1231        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1232      }
1233      return &m_samples;
1234    }
1235    
1236    
1237    const DataTypes::ValueType*
1238  #define PROC_OP(TYPE,X)                               \  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1239      for (int i=0;i<steps;++i,resultp+=resultStep) \  {
1240      { \    if (m_readytype!='E')
1241         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    {
1242         lroffset+=leftStep; \      throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1243         rroffset+=rightStep; \    }
1244      if (m_op==IDENTITY)
1245      {
1246        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1247      }
1248      size_t subroffset;
1249      size_t offset;
1250      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1251      roffset=m_samplesize*tid;
1252      offset=roffset;
1253      size_t loop=0;
1254      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1255      size_t outstep=getNoValues();
1256      size_t instep=m_left->getNoValues();
1257      switch (m_op)
1258      {
1259        case SWAP:
1260        for (loop=0;loop<numsteps;++loop)
1261        {
1262                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1263            subroffset+=instep;
1264            offset+=outstep;
1265      }      }
1266        break;
1267        default:
1268        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1269      }
1270      return &m_samples;
1271    }
1272    
1273    
1274    
 /*  
   \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.  
 */  
1275  // 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
1276  // 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.
1277  // 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 842  DataLazy::resolveNP1OUT(ValueType& v, si Line 1281  DataLazy::resolveNP1OUT(ValueType& v, si
1281  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1282  // For example, 2+Vector.  // For example, 2+Vector.
1283  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1284  DataTypes::ValueType*  const DataTypes::ValueType*
1285  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1286  {  {
1287  cout << "Resolve binary: " << toString() << endl;  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1288    
1289    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
1290      // first work out which of the children are expanded      // first work out which of the children are expanded
1291    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1292    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1293    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if (!leftExp && !rightExp)
1294    int steps=(bigloops?1:getNumDPPSample());    {
1295    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'.");
1296    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1297    {    bool leftScalar=(m_left->getRank()==0);
1298      EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");    bool rightScalar=(m_right->getRank()==0);
1299      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1300      chunksize=1;    // for scalar    {
1301    }          throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1302    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    }
1303    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    size_t leftsize=m_left->getNoValues();
1304    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    size_t rightsize=m_right->getNoValues();
1305      size_t chunksize=1;           // how many doubles will be processed in one go
1306      int leftstep=0;       // how far should the left offset advance after each step
1307      int rightstep=0;
1308      int numsteps=0;       // total number of steps for the inner loop
1309      int oleftstep=0;  // the o variables refer to the outer loop
1310      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1311      int onumsteps=1;
1312      
1313      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1314      bool RES=(rightExp && rightScalar);
1315      bool LS=(!leftExp && leftScalar); // left is a single scalar
1316      bool RS=(!rightExp && rightScalar);
1317      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1318      bool RN=(!rightExp && !rightScalar);
1319      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1320      bool REN=(rightExp && !rightScalar);
1321    
1322      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1323      {
1324        chunksize=m_left->getNumDPPSample()*leftsize;
1325        leftstep=0;
1326        rightstep=0;
1327        numsteps=1;
1328      }
1329      else if (LES || RES)
1330      {
1331        chunksize=1;
1332        if (LES)        // left is an expanded scalar
1333        {
1334            if (RS)
1335            {
1336               leftstep=1;
1337               rightstep=0;
1338               numsteps=m_left->getNumDPPSample();
1339            }
1340            else        // RN or REN
1341            {
1342               leftstep=0;
1343               oleftstep=1;
1344               rightstep=1;
1345               orightstep=(RN ? -(int)rightsize : 0);
1346               numsteps=rightsize;
1347               onumsteps=m_left->getNumDPPSample();
1348            }
1349        }
1350        else        // right is an expanded scalar
1351        {
1352            if (LS)
1353            {
1354               rightstep=1;
1355               leftstep=0;
1356               numsteps=m_right->getNumDPPSample();
1357            }
1358            else
1359            {
1360               rightstep=0;
1361               orightstep=1;
1362               leftstep=1;
1363               oleftstep=(LN ? -(int)leftsize : 0);
1364               numsteps=leftsize;
1365               onumsteps=m_right->getNumDPPSample();
1366            }
1367        }
1368      }
1369      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1370      {
1371        if (LEN)    // and Right will be a single value
1372        {
1373            chunksize=rightsize;
1374            leftstep=rightsize;
1375            rightstep=0;
1376            numsteps=m_left->getNumDPPSample();
1377            if (RS)
1378            {
1379               numsteps*=leftsize;
1380            }
1381        }
1382        else    // REN
1383        {
1384            chunksize=leftsize;
1385            rightstep=leftsize;
1386            leftstep=0;
1387            numsteps=m_right->getNumDPPSample();
1388            if (LS)
1389            {
1390               numsteps*=rightsize;
1391            }
1392        }
1393      }
1394    
1395      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1396      // Get the values of sub-expressions      // Get the values of sub-expressions
1397    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
1398    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1399      // the right child starts further along.  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1400    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;)
1401    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1402    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1403    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1404    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1405    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1406    
1407    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1408    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1409    
1410    
1411      roffset=m_samplesize*tid;
1412      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
1413    switch(m_op)    switch(m_op)
1414    {    {
1415      case ADD:      case ADD:
# Line 888  cout << "Resolve binary: " << toString() Line 1430  cout << "Resolve binary: " << toString()
1430      default:      default:
1431      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1432    }    }
1433    roffset=offset;    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1434    return &v;    return &m_samples;
1435  }  }
1436    
1437    
 /*  
   \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.  
 */  
1438  // 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
1439  // 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.
1440  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1441  DataTypes::ValueType*  const DataTypes::ValueType*
1442  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
1443  {  {
1444  cout << "Resolve TensorProduct: " << toString() << endl;  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1445    
1446    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
1447      // first work out which of the children are expanded      // first work out which of the children are expanded
1448    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1449    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1450    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1451    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1452    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1453    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0  
1454      // Get the values of sub-expressions (leave a gap of one sample for the result).    int resultStep=getNoValues();
1455    const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);    roffset=m_samplesize*tid;
1456    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);    size_t offset=roffset;
1457    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved  
1458      const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1459    
1460      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1461    
1462    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1463    cout << getNoValues() << endl;)
1464    
1465    
1466    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1467    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1468    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1469    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1470    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1471    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1472    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1473    
1474      double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
1475    switch(m_op)    switch(m_op)
1476    {    {
1477      case PROD:      case PROD:
# Line 932  cout << "Resolve TensorProduct: " << toS Line 1479  cout << "Resolve TensorProduct: " << toS
1479      {      {
1480            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1481            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1482    
1483    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1484    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1485    
1486            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);
1487    
1488        lroffset+=leftStep;        lroffset+=leftStep;
1489        rroffset+=rightStep;        rroffset+=rightStep;
1490      }      }
# Line 941  cout << "Resolve TensorProduct: " << toS Line 1493  cout << "Resolve TensorProduct: " << toS
1493      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1494    }    }
1495    roffset=offset;    roffset=offset;
1496    return &v;    return &m_samples;
1497  }  }
1498    
1499    
1500    const DataTypes::ValueType*
1501    DataLazy::resolveSample(int sampleNo, size_t& roffset)
1502    {
1503    #ifdef _OPENMP
1504        int tid=omp_get_thread_num();
1505    #else
1506        int tid=0;
1507    #endif
1508        return resolveNodeSample(tid, sampleNo, roffset);
1509    }
1510    
 /*  
   \brief Compute the value of the expression 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.  
1511    
1512    The return value will be an existing vector so do not deallocate it.  // This needs to do the work of the identity constructor
1513  */  void
1514  // the vector and the offset are a place where the method could write its data if it wishes  DataLazy::resolveToIdentity()
1515  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  {
1516  // Hence the return value to indicate where the data is actually stored.     if (m_op==IDENTITY)
1517  // Regardless, the storage should be assumed to be used, even if it isn't.      return;
1518       DataReady_ptr p=resolveNodeWorker();
1519       makeIdentity(p);
1520    }
1521    
1522  // the roffset is the offset within the returned vector where the data begins  void DataLazy::makeIdentity(const DataReady_ptr& p)
 const DataTypes::ValueType*  
 DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  
1523  {  {
1524  cout << "Resolve sample " << toString() << endl;     m_op=IDENTITY;
1525      // collapse so we have a 'E' node or an IDENTITY for some other type     m_axis_offset=0;
1526    if (m_readytype!='E' && m_op!=IDENTITY)     m_transpose=0;
1527    {     m_SL=m_SM=m_SR=0;
1528      collapse();     m_children=m_height=0;
1529    }     m_id=p;
1530    if (m_op==IDENTITY)       if(p->isConstant()) {m_readytype='C';}
1531    {     else if(p->isExpanded()) {m_readytype='E';}
1532      const ValueType& vec=m_id->getVector();     else if (p->isTagged()) {m_readytype='T';}
1533      if (m_readytype=='C')     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1534      {     m_samplesize=p->getNumDPPSample()*p->getNoValues();
1535      roffset=0;     m_left.reset();
1536      return &(vec);     m_right.reset();
     }  
     roffset=m_id->getPointOffset(sampleNo, 0);  
     return &(vec);  
   }  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer Error - Collapse did not produce an expanded node.");  
   }  
   switch (getOpgroup(m_op))  
   {  
   case G_UNARY: 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_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);  
   default:  
     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
   }  
1537  }  }
1538    
1539    
 // 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.  
1540  DataReady_ptr  DataReady_ptr
1541  DataLazy::resolve()  DataLazy::resolve()
1542  {  {
1543        resolveToIdentity();
1544        return m_id;
1545    }
1546    
1547  cout << "Sample size=" << m_samplesize << endl;  // This version of resolve uses storage in each node to hold results
1548  cout << "Buffers=" << m_buffsRequired << endl;  DataReady_ptr
1549    DataLazy::resolveNodeWorker()
1550    {
1551    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
1552    {    {
1553      collapse();      collapse();
# Line 1016  cout << "Buffers=" << m_buffsRequired << Line 1557  cout << "Buffers=" << m_buffsRequired <<
1557      return m_id;      return m_id;
1558    }    }
1559      // 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);  
 cout << "Buffer created with size=" << v.size() << endl;  
1560    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1561    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
1562    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1563    
1564    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1565    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1566    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Storage for answer
1567    size_t resoffset=0;       // where in the vector to find the answer  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1568    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,res) schedule(static)
1569    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1570    {    {
1571  cout << "################################# " << sample << endl;      size_t roffset=0;
1572  #ifdef _OPENMP  #ifdef _OPENMP
1573      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1574  #else  #else
1575      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.      res=resolveNodeSample(0,sample,roffset);
1576  #endif  #endif
1577  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1578      outoffset=result->getPointOffset(sample,0);  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1579  cerr << "offset=" << outoffset << endl;      DataVector::size_type outoffset=result->getPointOffset(sample,0);
1580      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));
     {  
     resvec[outoffset]=(*res)[resoffset];  
     }  
 cerr << "*********************************" << endl;  
1581    }    }
1582    return resptr;    return resptr;
1583  }  }
# Line 1058  DataLazy::toString() const Line 1587  DataLazy::toString() const
1587  {  {
1588    ostringstream oss;    ostringstream oss;
1589    oss << "Lazy Data:";    oss << "Lazy Data:";
1590    intoString(oss);    if (escriptParams.getPRINT_LAZY_TREE()==0)
1591      {
1592          intoString(oss);
1593      }
1594      else
1595      {
1596        oss << endl;
1597        intoTreeString(oss,"");
1598      }
1599    return oss.str();    return oss.str();
1600  }  }
1601    
# Line 1066  DataLazy::toString() const Line 1603  DataLazy::toString() const
1603  void  void
1604  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1605  {  {
1606    //    oss << "[" << m_children <<";"<<m_height <<"]";
1607    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1608    {    {
1609    case G_IDENTITY:    case G_IDENTITY:
# Line 1095  DataLazy::intoString(ostringstream& oss) Line 1633  DataLazy::intoString(ostringstream& oss)
1633      oss << ')';      oss << ')';
1634      break;      break;
1635    case G_UNARY:    case G_UNARY:
1636      case G_UNARY_P:
1637    case G_NP1OUT:    case G_NP1OUT:
1638      case G_NP1OUT_P:
1639      case G_REDUCTION:
1640      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1641      m_left->intoString(oss);      m_left->intoString(oss);
1642      oss << ')';      oss << ')';
# Line 1107  DataLazy::intoString(ostringstream& oss) Line 1648  DataLazy::intoString(ostringstream& oss)
1648      m_right->intoString(oss);      m_right->intoString(oss);
1649      oss << ')';      oss << ')';
1650      break;      break;
1651      case G_NP1OUT_2P:
1652        oss << opToString(m_op) << '(';
1653        m_left->intoString(oss);
1654        oss << ", " << m_axis_offset << ", " << m_transpose;
1655        oss << ')';
1656        break;
1657    default:    default:
1658      oss << "UNKNOWN";      oss << "UNKNOWN";
1659    }    }
1660  }  }
1661    
1662    
1663    void
1664    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1665    {
1666      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1667      switch (getOpgroup(m_op))
1668      {
1669      case G_IDENTITY:
1670        if (m_id->isExpanded())
1671        {
1672           oss << "E";
1673        }
1674        else if (m_id->isTagged())
1675        {
1676          oss << "T";
1677        }
1678        else if (m_id->isConstant())
1679        {
1680          oss << "C";
1681        }
1682        else
1683        {
1684          oss << "?";
1685        }
1686        oss << '@' << m_id.get() << endl;
1687        break;
1688      case G_BINARY:
1689        oss << opToString(m_op) << endl;
1690        indent+='.';
1691        m_left->intoTreeString(oss, indent);
1692        m_right->intoTreeString(oss, indent);
1693        break;
1694      case G_UNARY:
1695      case G_UNARY_P:
1696      case G_NP1OUT:
1697      case G_NP1OUT_P:
1698      case G_REDUCTION:
1699        oss << opToString(m_op) << endl;
1700        indent+='.';
1701        m_left->intoTreeString(oss, indent);
1702        break;
1703      case G_TENSORPROD:
1704        oss << opToString(m_op) << endl;
1705        indent+='.';
1706        m_left->intoTreeString(oss, indent);
1707        m_right->intoTreeString(oss, indent);
1708        break;
1709      case G_NP1OUT_2P:
1710        oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1711        indent+='.';
1712        m_left->intoTreeString(oss, indent);
1713        break;
1714      default:
1715        oss << "UNKNOWN";
1716      }
1717    }
1718    
1719    
1720  DataAbstract*  DataAbstract*
1721  DataLazy::deepCopy()  DataLazy::deepCopy()
1722  {  {
1723    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1724    {    {
1725    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1726    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY:
1727      case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1728      case G_UNARY_P:   return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1729    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);
1730    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);
1731    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);
1732      case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
1733      case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1734    default:    default:
1735      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)+".");
1736    }    }
1737  }  }
1738    
1739    
1740    
1741  // There is no single, natural interpretation of getLength on DataLazy.  // There is no single, natural interpretation of getLength on DataLazy.
1742  // Instances of DataReady can look at the size of their vectors.  // Instances of DataReady can look at the size of their vectors.
1743  // 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 1203  DataLazy::getPointOffset(int sampleNo, Line 1813  DataLazy::getPointOffset(int sampleNo,
1813    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).");
1814  }  }
1815    
1816  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
1817  // 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.  
1818  void  void
1819  DataLazy::setToZero()  DataLazy::setToZero()
1820  {  {
1821    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
1822    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1823    m_op=IDENTITY;  //   m_op=IDENTITY;
1824    m_right.reset();    //   m_right.reset();  
1825    m_left.reset();  //   m_left.reset();
1826    m_readytype='C';  //   m_readytype='C';
1827    m_buffsRequired=1;  //   m_buffsRequired=1;
1828    
1829      privdebug=privdebug;  // to stop the compiler complaining about unused privdebug
1830      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1831    }
1832    
1833    bool
1834    DataLazy::actsExpanded() const
1835    {
1836        return (m_readytype=='E');
1837  }  }
1838    
1839  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26