/[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 2166 by jfenwick, Tue Dec 16 06:08:02 2008 UTC revision 2514 by jfenwick, Fri Jul 3 00:57:45 2009 UTC
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31    #include "EscriptParams.h"
32    
33    #include <iomanip>      // for some fancy formatting in debug
34    
35  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
36  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
37  namespace  namespace
# Line 38  bool privdebug=false; Line 42  bool privdebug=false;
42  #define DISABLEDEBUG privdebug=false;  #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 109  enum ES_opgroup Line 118  enum ES_opgroup
118     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     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_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
121     G_TENSORPROD     // general tensor product     G_TENSORPROD,    // general tensor product
122       G_NP1OUT_2P      // non-pointwise op with one output requiring two params
123  };  };
124    
125    
# Line 123  string ES_opstrings[]={"UNKNOWN","IDENTI Line 133  string ES_opstrings[]={"UNKNOWN","IDENTI
133              "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",
134              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
135              "prod",              "prod",
136              "transpose", "trace"};              "transpose", "trace",
137  int ES_opcount=40;              "swapaxes"};
138    int ES_opcount=41;
139  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,
140              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
141              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 133  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT Line 144  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT
144              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
145              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
146              G_TENSORPROD,              G_TENSORPROD,
147              G_NP1OUT_P, G_NP1OUT_P};              G_NP1OUT_P, G_NP1OUT_P,
148                G_NP1OUT_2P};
149  inline  inline
150  ES_opgroup  ES_opgroup
151  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 194  resultShape(DataAbstract_ptr left, DataA Line 206  resultShape(DataAbstract_ptr left, DataA
206  // return the shape for "op left"  // return the shape for "op left"
207    
208  DataTypes::ShapeType  DataTypes::ShapeType
209  resultShape(DataAbstract_ptr left, ES_optype op)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
210  {  {
211      switch(op)      switch(op)
212      {      {
213          case TRANS:          case TRANS:
214          return left->getShape();         {            // for the scoping of variables
215            const DataTypes::ShapeType& s=left->getShape();
216            DataTypes::ShapeType sh;
217            int rank=left->getRank();
218            if (axis_offset<0 || axis_offset>rank)
219            {
220                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
221                }
222                for (int i=0; i<rank; i++)
223            {
224               int index = (axis_offset+i)%rank;
225                   sh.push_back(s[index]); // Append to new shape
226                }
227            return sh;
228           }
229      break;      break;
230      case TRACE:      case TRACE:
231          return DataTypes::scalarShape;         {
232            int rank=left->getRank();
233            if (rank<2)
234            {
235               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
236            }
237            if ((axis_offset>rank-2) || (axis_offset<0))
238            {
239               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
240            }
241            if (rank==2)
242            {
243               return DataTypes::scalarShape;
244            }
245            else if (rank==3)
246            {
247               DataTypes::ShapeType sh;
248                   if (axis_offset==0)
249               {
250                    sh.push_back(left->getShape()[2]);
251                   }
252                   else     // offset==1
253               {
254                sh.push_back(left->getShape()[0]);
255                   }
256               return sh;
257            }
258            else if (rank==4)
259            {
260               DataTypes::ShapeType sh;
261               const DataTypes::ShapeType& s=left->getShape();
262                   if (axis_offset==0)
263               {
264                    sh.push_back(s[2]);
265                    sh.push_back(s[3]);
266                   }
267                   else if (axis_offset==1)
268               {
269                    sh.push_back(s[0]);
270                    sh.push_back(s[3]);
271                   }
272               else     // offset==2
273               {
274                sh.push_back(s[0]);
275                sh.push_back(s[1]);
276               }
277               return sh;
278            }
279            else        // unknown rank
280            {
281               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
282            }
283           }
284      break;      break;
285          default:          default:
286      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)+".");
287      }      }
288  }  }
289    
290    DataTypes::ShapeType
291    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
292    {
293         // This code taken from the Data.cpp swapaxes() method
294         // Some of the checks are probably redundant here
295         int axis0_tmp,axis1_tmp;
296         const DataTypes::ShapeType& s=left->getShape();
297         DataTypes::ShapeType out_shape;
298         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
299         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
300         int rank=left->getRank();
301         if (rank<2) {
302            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
303         }
304         if (axis0<0 || axis0>rank-1) {
305            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
306         }
307         if (axis1<0 || axis1>rank-1) {
308             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
309         }
310         if (axis0 == axis1) {
311             throw DataException("Error - Data::swapaxes: axis indices must be different.");
312         }
313         if (axis0 > axis1) {
314             axis0_tmp=axis1;
315             axis1_tmp=axis0;
316         } else {
317             axis0_tmp=axis0;
318             axis1_tmp=axis1;
319         }
320         for (int i=0; i<rank; i++) {
321           if (i == axis0_tmp) {
322              out_shape.push_back(s[axis1_tmp]);
323           } else if (i == axis1_tmp) {
324              out_shape.push_back(s[axis0_tmp]);
325           } else {
326              out_shape.push_back(s[i]);
327           }
328         }
329        return out_shape;
330    }
331    
332    
333  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
334  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
335  // 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 273  GTPShape(DataAbstract_ptr left, DataAbst Line 394  GTPShape(DataAbstract_ptr left, DataAbst
394    return shape2;    return shape2;
395  }  }
396    
   
 // 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)+".");  
 //    }  
 // }  
   
397  // determine the number of samples requires to evaluate an expression combining left and right  // determine the number of samples requires to evaluate an expression combining left and right
398  // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.  // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
399  // The same goes for G_TENSORPROD  // The same goes for G_TENSORPROD
400  // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefths.  // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
401  // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined  // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
402  // multiple times  // multiple times
403  int  int
# Line 307  calcBuffs(const DataLazy_ptr& left, cons Line 412  calcBuffs(const DataLazy_ptr& left, cons
412     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
413     case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);     case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
414     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
415       case G_NP1OUT_2P: return 1+max(left->getBuffsRequired(),1);
416     default:     default:
417      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
418     }     }
# Line 328  opToString(ES_optype op) Line 434  opToString(ES_optype op)
434    return ES_opstrings[op];    return ES_opstrings[op];
435  }  }
436    
437    #ifdef LAZY_NODE_STORAGE
438    void DataLazy::LazyNodeSetup()
439    {
440    #ifdef _OPENMP
441        int numthreads=omp_get_max_threads();
442        m_samples.resize(numthreads*m_samplesize);
443        m_sampleids=new int[numthreads];
444        for (int i=0;i<numthreads;++i)
445        {
446            m_sampleids[i]=-1;  
447        }
448    #else
449        m_samples.resize(m_samplesize);
450        m_sampleids=new int[1];
451        m_sampleids[0]=-1;
452    #endif  // _OPENMP
453    }
454    #endif   // LAZY_NODE_STORAGE
455    
456    
457    // Creates an identity node
458  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
459      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
460      m_op(IDENTITY),  #ifdef LAZY_NODE_STORAGE
461      m_axis_offset(0),      ,m_sampleids(0),
462      m_transpose(0),      m_samples(1)
463      m_SL(0), m_SM(0), m_SR(0)  #endif
464  {  {
465     if (p->isLazy())     if (p->isLazy())
466     {     {
# Line 345  DataLazy::DataLazy(DataAbstract_ptr p) Line 471  DataLazy::DataLazy(DataAbstract_ptr p)
471     }     }
472     else     else
473     {     {
474      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
475      if(p->isConstant()) {m_readytype='C';}      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
476      else if(p->isExpanded()) {m_readytype='E';}      makeIdentity(dr);
477      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.");}  
478     }     }
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_maxsamplesize=m_samplesize;  
479  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
480  }  }
481    
   
   
   
482  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
483      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
484      m_op(op),      m_op(op),
# Line 386  DataLazy::DataLazy(DataAbstract_ptr left Line 505  DataLazy::DataLazy(DataAbstract_ptr left
505     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
506     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
507     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
508       m_children=m_left->m_children+1;
509       m_height=m_left->m_height+1;
510    #ifdef LAZY_NODE_STORAGE
511       LazyNodeSetup();
512    #endif
513       SIZELIMIT
514  }  }
515    
516    
# Line 395  DataLazy::DataLazy(DataAbstract_ptr left Line 520  DataLazy::DataLazy(DataAbstract_ptr left
520      m_op(op),      m_op(op),
521      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
522  {  {
523    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
524     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
525     {     {
526      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 411  DataLazy::DataLazy(DataAbstract_ptr left Line 537  DataLazy::DataLazy(DataAbstract_ptr left
537     {     {
538      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
539      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
540    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
541     }     }
542     left->operandCheck(*right);     left->operandCheck(*right);
543    
544     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
545     {     {
546      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
547    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
548     }     }
549     else     else
550     {     {
551      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
552    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
553     }     }
554     if (right->isLazy())     if (right->isLazy())
555     {     {
556      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
557    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
558     }     }
559     else     else
560     {     {
561      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
562    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
563     }     }
564     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
565     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 447  DataLazy::DataLazy(DataAbstract_ptr left Line 578  DataLazy::DataLazy(DataAbstract_ptr left
578     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
579     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
580     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
581       m_children=m_left->m_children+m_right->m_children+2;
582       m_height=max(m_left->m_height,m_right->m_height)+1;
583    #ifdef LAZY_NODE_STORAGE
584       LazyNodeSetup();
585    #endif
586       SIZELIMIT
587  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
588  }  }
589    
# Line 476  DataLazy::DataLazy(DataAbstract_ptr left Line 613  DataLazy::DataLazy(DataAbstract_ptr left
613      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
614      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
615     }     }
616     left->operandCheck(*right);  //    left->operandCheck(*right);
617    
618     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
619     {     {
# Line 511  DataLazy::DataLazy(DataAbstract_ptr left Line 648  DataLazy::DataLazy(DataAbstract_ptr left
648     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
649     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
650     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
651       m_children=m_left->m_children+m_right->m_children+2;
652       m_height=max(m_left->m_height,m_right->m_height)+1;
653    #ifdef LAZY_NODE_STORAGE
654       LazyNodeSetup();
655    #endif
656       SIZELIMIT
657  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
658  }  }
659    
660    
661  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
662      : parent(left->getFunctionSpace(), resultShape(left,op)),      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
663      m_op(op),      m_op(op),
664      m_axis_offset(axis_offset),      m_axis_offset(axis_offset),
665      m_transpose(0),      m_transpose(0),
# Line 540  DataLazy::DataLazy(DataAbstract_ptr left Line 683  DataLazy::DataLazy(DataAbstract_ptr left
683     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
684     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
685     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
686       m_children=m_left->m_children+1;
687       m_height=m_left->m_height+1;
688    #ifdef LAZY_NODE_STORAGE
689       LazyNodeSetup();
690    #endif
691       SIZELIMIT
692  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
693  }  }
694    
# Line 568  DataLazy::DataLazy(DataAbstract_ptr left Line 717  DataLazy::DataLazy(DataAbstract_ptr left
717     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
718     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
719     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
720       m_children=m_left->m_children+1;
721       m_height=m_left->m_height+1;
722    #ifdef LAZY_NODE_STORAGE
723       LazyNodeSetup();
724    #endif
725       SIZELIMIT
726  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
727  }  }
728    
729    
730    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
731        : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
732        m_op(op),
733        m_axis_offset(axis0),
734        m_transpose(axis1),
735        m_tol(0)
736    {
737       if ((getOpgroup(op)!=G_NP1OUT_2P))
738       {
739        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
740       }
741       DataLazy_ptr lleft;
742       if (!left->isLazy())
743       {
744        lleft=DataLazy_ptr(new DataLazy(left));
745       }
746       else
747       {
748        lleft=dynamic_pointer_cast<DataLazy>(left);
749       }
750       m_readytype=lleft->m_readytype;
751       m_left=lleft;
752       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
753       m_samplesize=getNumDPPSample()*getNoValues();
754       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
755       m_children=m_left->m_children+1;
756       m_height=m_left->m_height+1;
757    #ifdef LAZY_NODE_STORAGE
758       LazyNodeSetup();
759    #endif
760       SIZELIMIT
761    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
762    }
763    
764  DataLazy::~DataLazy()  DataLazy::~DataLazy()
765  {  {
766    #ifdef LAZY_NODE_SETUP
767       delete[] m_sampleids;
768       delete[] m_samples;
769    #endif
770  }  }
771    
772    
# Line 589  DataLazy::getMaxSampleSize() const Line 783  DataLazy::getMaxSampleSize() const
783      return m_maxsamplesize;      return m_maxsamplesize;
784  }  }
785    
786    
787    
788    size_t
789    DataLazy::getSampleBufferSize() const
790    {
791        return m_maxsamplesize*(max(1,m_buffsRequired));
792    }
793    
794  /*  /*
795    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
796    This does the work for the collapse method.    This does the work for the collapse method.
# Line 728  DataLazy::collapseToReady() Line 930  DataLazy::collapseToReady()
930      case TRACE:      case TRACE:
931      result=left.trace(m_axis_offset);      result=left.trace(m_axis_offset);
932      break;      break;
933        case SWAP:
934        result=left.swapaxes(m_axis_offset, m_transpose);
935        break;
936      default:      default:
937      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)+".");
938    }    }
# Line 1023  LAZYDEBUG(cerr << "Result of trace=" << Line 1228  LAZYDEBUG(cerr << "Result of trace=" <<
1228  }  }
1229    
1230    
1231    /*
1232      \brief Compute the value of the expression (unary operation with int params) for the given sample.
1233      \return Vector which stores the value of the subexpression for the given sample.
1234      \param v A vector to store intermediate results.
1235      \param offset Index in v to begin storing results.
1236      \param sampleNo Sample number to evaluate.
1237      \param roffset (output parameter) the offset in the return vector where the result begins.
1238    
1239      The return value will be an existing vector so do not deallocate it.
1240      If the result is stored in v it should be stored at the offset given.
1241      Everything from offset to the end of v should be considered available for this method to use.
1242    */
1243    DataTypes::ValueType*
1244    DataLazy::resolveNP1OUT_2P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1245    {
1246        // we assume that any collapsing has been done before we get here
1247        // since we only have one argument we don't need to think about only
1248        // processing single points.
1249      if (m_readytype!='E')
1250      {
1251        throw DataException("Programmer error - resolveNP1OUT_2P should only be called on expanded Data.");
1252      }
1253        // since we can't write the result over the input, we need a result offset further along
1254      size_t subroffset;
1255      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1256    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1257    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1258      roffset=offset;
1259      size_t loop=0;
1260      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1261      size_t outstep=getNoValues();
1262      size_t instep=m_left->getNoValues();
1263    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1264      switch (m_op)
1265      {
1266        case SWAP:
1267        for (loop=0;loop<numsteps;++loop)
1268        {
1269                DataMaths::swapaxes(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset, m_transpose);
1270            subroffset+=instep;
1271            offset+=outstep;
1272        }
1273        break;
1274        default:
1275        throw DataException("Programmer error - resolveNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1276      }
1277      return &v;
1278    }
1279    
1280    
1281    
1282  #define PROC_OP(TYPE,X)                               \  #define PROC_OP(TYPE,X)                               \
1283      for (int j=0;j<onumsteps;++j)\      for (int j=0;j<onumsteps;++j)\
1284      {\      {\
# Line 1121  LAZYDEBUG(cout << "Resolve binary: " << Line 1377  LAZYDEBUG(cout << "Resolve binary: " <<
1377             leftstep=0;             leftstep=0;
1378             oleftstep=1;             oleftstep=1;
1379             rightstep=1;             rightstep=1;
1380             orightstep=(RN?-rightsize:0);             orightstep=(RN ? -(int)rightsize : 0);
1381             numsteps=rightsize;             numsteps=rightsize;
1382             onumsteps=m_left->getNumDPPSample();             onumsteps=m_left->getNumDPPSample();
1383          }          }
# Line 1139  LAZYDEBUG(cout << "Resolve binary: " << Line 1395  LAZYDEBUG(cout << "Resolve binary: " <<
1395             rightstep=0;             rightstep=0;
1396             orightstep=1;             orightstep=1;
1397             leftstep=1;             leftstep=1;
1398             oleftstep=(LN?-leftsize:0);             oleftstep=(LN ? -(int)leftsize : 0);
1399             numsteps=leftsize;             numsteps=leftsize;
1400             onumsteps=m_right->getNumDPPSample();             onumsteps=m_right->getNumDPPSample();
1401          }          }
# Line 1184  LAZYDEBUG(cout << " numsteps=" << numste Line 1440  LAZYDEBUG(cout << " numsteps=" << numste
1440  LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)  LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1441  LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)  LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1442  LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)  LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1443    
1444    
1445    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1446    switch(m_op)    switch(m_op)
1447    {    {
# Line 1205  LAZYDEBUG(cout << "" << LS << RS << LN < Line 1463  LAZYDEBUG(cout << "" << LS << RS << LN <
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;      roffset=offset;
1467    return &v;    return &v;
1468  }  }
1469    
# Line 1229  LAZYDEBUG(cout << "" << LS << RS << LN < Line 1487  LAZYDEBUG(cout << "" << LS << RS << LN <
1487  DataTypes::ValueType*  DataTypes::ValueType*
1488  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1489  {  {
1490  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1491    
1492    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
1493      // first work out which of the children are expanded      // first work out which of the children are expanded
1494    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1495    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1496    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1497    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);  /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1498    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1499    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1500      int rightStep=(rightExp?m_right->getNoValues() : 0);
1501    
1502      int resultStep=getNoValues();
1503      // Get the values of sub-expressions (leave a gap of one sample for the result).      // Get the values of sub-expressions (leave a gap of one sample for the result).
1504    int gap=offset+m_left->getMaxSampleSize();    // actually only needs to be m_left->m_samplesize    int gap=offset+m_samplesize;  
1505    
1506    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1507    
1508    const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1509    gap+=m_right->getMaxSampleSize();    gap+=m_left->getMaxSampleSize();
1510    const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);  
1511  LAZYDEBUG(cout << "Post sub calls: " << toString() << endl;)  
1512    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1513    
1514    
1515      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1516    
1517    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1518    cout << getNoValues() << endl;)
1519    LAZYDEBUG(cerr << "Result of left=";)
1520    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1521    
1522    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1523    {
1524    cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1525    if (i%4==0) cout << endl;
1526    })
1527    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1528    LAZYDEBUG(
1529    for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1530    {
1531    cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1532    if (i%4==0) cout << endl;
1533    }
1534    cerr << endl;
1535    )
1536    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1537  LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)  LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1538  LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)  LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1539  LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)  LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1540  LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)  LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1541  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1542    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1543    
1544    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1545    switch(m_op)    switch(m_op)
1546    {    {
1547      case PROD:      case PROD:
1548      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1549      {      {
1550    
1551  LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)  LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1552  LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)  LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1553  LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)  LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1554    
1555            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1556            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1557    
1558  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1559  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1560    
1561            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);
1562    
1563    LAZYDEBUG(cout << "Results--\n";
1564    {
1565      DataVector dv(getNoValues());
1566    for (int z=0;z<getNoValues();++z)
1567    {
1568      cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1569      if (z%4==0) cout << endl;
1570      dv[z]=resultp[z];
1571    }
1572    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1573    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1574    }
1575    )
1576        lroffset+=leftStep;        lroffset+=leftStep;
1577        rroffset+=rightStep;        rroffset+=rightStep;
1578      }      }
# Line 1276  LAZYDEBUG(cout << DataTypes::pointToStri Line 1585  LAZYDEBUG(cout << DataTypes::pointToStri
1585  }  }
1586    
1587    
1588    #ifdef LAZY_NODE_STORAGE
1589    
1590    // The result will be stored in m_samples
1591    // The return value is a pointer to the DataVector, offset is the offset within the return value
1592    const DataTypes::ValueType*
1593    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
1594    {
1595    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1596        // collapse so we have a 'E' node or an IDENTITY for some other type
1597      if (m_readytype!='E' && m_op!=IDENTITY)
1598      {
1599        collapse();
1600      }
1601      if (m_op==IDENTITY)  
1602      {
1603        const ValueType& vec=m_id->getVectorRO();
1604    //     if (m_readytype=='C')
1605    //     {
1606    //  roffset=0;      // all samples read from the same position
1607    //  return &(m_samples);
1608    //     }
1609        roffset=m_id->getPointOffset(sampleNo, 0);
1610        return &(vec);
1611      }
1612      if (m_readytype!='E')
1613      {
1614        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1615      }
1616      if (m_sampleids[tid]==sampleNo)
1617      {
1618        roffset=tid*m_samplesize;
1619        return &(m_samples);        // sample is already resolved
1620      }
1621      m_sampleids[tid]=sampleNo;
1622      switch (getOpgroup(m_op))
1623      {
1624      case G_UNARY:
1625      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1626      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1627      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1628      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1629      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1630      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1631      default:
1632        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1633      }
1634    }
1635    
1636    const DataTypes::ValueType*
1637    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
1638    {
1639        // we assume that any collapsing has been done before we get here
1640        // since we only have one argument we don't need to think about only
1641        // processing single points.
1642        // we will also know we won't get identity nodes
1643      if (m_readytype!='E')
1644      {
1645        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1646      }
1647      if (m_op==IDENTITY)
1648      {
1649        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1650      }
1651      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1652      const double* left=&((*leftres)[roffset]);
1653      roffset=m_samplesize*tid;
1654      double* result=&(m_samples[roffset]);
1655      switch (m_op)
1656      {
1657        case SIN:  
1658        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
1659        break;
1660        case COS:
1661        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
1662        break;
1663        case TAN:
1664        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
1665        break;
1666        case ASIN:
1667        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1668        break;
1669        case ACOS:
1670        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1671        break;
1672        case ATAN:
1673        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1674        break;
1675        case SINH:
1676        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1677        break;
1678        case COSH:
1679        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1680        break;
1681        case TANH:
1682        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1683        break;
1684        case ERF:
1685    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1686        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1687    #else
1688        tensor_unary_operation(m_samplesize, left, result, ::erf);
1689        break;
1690    #endif
1691       case ASINH:
1692    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1693        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1694    #else
1695        tensor_unary_operation(m_samplesize, left, result, ::asinh);
1696    #endif  
1697        break;
1698       case ACOSH:
1699    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1700        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1701    #else
1702        tensor_unary_operation(m_samplesize, left, result, ::acosh);
1703    #endif  
1704        break;
1705       case ATANH:
1706    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1707        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1708    #else
1709        tensor_unary_operation(m_samplesize, left, result, ::atanh);
1710    #endif  
1711        break;
1712        case LOG10:
1713        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1714        break;
1715        case LOG:
1716        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1717        break;
1718        case SIGN:
1719        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1720        break;
1721        case ABS:
1722        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1723        break;
1724        case NEG:
1725        tensor_unary_operation(m_samplesize, left, result, negate<double>());
1726        break;
1727        case POS:
1728        // it doesn't mean anything for delayed.
1729        // it will just trigger a deep copy of the lazy object
1730        throw DataException("Programmer error - POS not supported for lazy data.");
1731        break;
1732        case EXP:
1733        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1734        break;
1735        case SQRT:
1736        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1737        break;
1738        case RECIP:
1739        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1740        break;
1741        case GZ:
1742        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1743        break;
1744        case LZ:
1745        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1746        break;
1747        case GEZ:
1748        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1749        break;
1750        case LEZ:
1751        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1752        break;
1753    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1754        case NEZ:
1755        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1756        break;
1757        case EZ:
1758        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1759        break;
1760    
1761        default:
1762        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1763      }
1764      return &(m_samples);
1765    }
1766    
1767    
1768    const DataTypes::ValueType*
1769    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
1770    {
1771        // we assume that any collapsing has been done before we get here
1772        // since we only have one argument we don't need to think about only
1773        // processing single points.
1774      if (m_readytype!='E')
1775      {
1776        throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1777      }
1778      if (m_op==IDENTITY)
1779      {
1780        throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1781      }
1782      size_t subroffset;
1783      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1784      roffset=m_samplesize*tid;
1785      size_t loop=0;
1786      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1787      size_t step=getNoValues();
1788      size_t offset=roffset;
1789      switch (m_op)
1790      {
1791        case SYM:
1792        for (loop=0;loop<numsteps;++loop)
1793        {
1794            DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1795            subroffset+=step;
1796            offset+=step;
1797        }
1798        break;
1799        case NSYM:
1800        for (loop=0;loop<numsteps;++loop)
1801        {
1802            DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1803            subroffset+=step;
1804            offset+=step;
1805        }
1806        break;
1807        default:
1808        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1809      }
1810      return &m_samples;
1811    }
1812    
1813    const DataTypes::ValueType*
1814    DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
1815    {
1816        // we assume that any collapsing has been done before we get here
1817        // since we only have one argument we don't need to think about only
1818        // processing single points.
1819      if (m_readytype!='E')
1820      {
1821        throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1822      }
1823      if (m_op==IDENTITY)
1824      {
1825        throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1826      }
1827      size_t subroffset;
1828      size_t offset;
1829      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1830      roffset=m_samplesize*tid;
1831      offset=roffset;
1832      size_t loop=0;
1833      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1834      size_t outstep=getNoValues();
1835      size_t instep=m_left->getNoValues();
1836      switch (m_op)
1837      {
1838        case TRACE:
1839        for (loop=0;loop<numsteps;++loop)
1840        {
1841                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1842            subroffset+=instep;
1843            offset+=outstep;
1844        }
1845        break;
1846        case TRANS:
1847        for (loop=0;loop<numsteps;++loop)
1848        {
1849                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1850            subroffset+=instep;
1851            offset+=outstep;
1852        }
1853        break;
1854        default:
1855        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1856      }
1857      return &m_samples;
1858    }
1859    
1860    
1861    const DataTypes::ValueType*
1862    DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1863    {
1864      if (m_readytype!='E')
1865      {
1866        throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1867      }
1868      if (m_op==IDENTITY)
1869      {
1870        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1871      }
1872      size_t subroffset;
1873      size_t offset;
1874      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1875      roffset=m_samplesize*tid;
1876      offset=roffset;
1877      size_t loop=0;
1878      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1879      size_t outstep=getNoValues();
1880      size_t instep=m_left->getNoValues();
1881      switch (m_op)
1882      {
1883        case SWAP:
1884        for (loop=0;loop<numsteps;++loop)
1885        {
1886                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1887            subroffset+=instep;
1888            offset+=outstep;
1889        }
1890        break;
1891        default:
1892        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1893      }
1894      return &m_samples;
1895    }
1896    
1897    
1898    
1899    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1900    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1901    // If both children are expanded, then we can process them in a single operation (we treat
1902    // the whole sample as one big datapoint.
1903    // If one of the children is not expanded, then we need to treat each point in the sample
1904    // individually.
1905    // There is an additional complication when scalar operations are considered.
1906    // For example, 2+Vector.
1907    // In this case each double within the point is treated individually
1908    const DataTypes::ValueType*
1909    DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1910    {
1911    LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1912    
1913      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1914        // first work out which of the children are expanded
1915      bool leftExp=(m_left->m_readytype=='E');
1916      bool rightExp=(m_right->m_readytype=='E');
1917      if (!leftExp && !rightExp)
1918      {
1919        throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1920      }
1921      bool leftScalar=(m_left->getRank()==0);
1922      bool rightScalar=(m_right->getRank()==0);
1923      if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1924      {
1925        throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1926      }
1927      size_t leftsize=m_left->getNoValues();
1928      size_t rightsize=m_right->getNoValues();
1929      size_t chunksize=1;           // how many doubles will be processed in one go
1930      int leftstep=0;       // how far should the left offset advance after each step
1931      int rightstep=0;
1932      int numsteps=0;       // total number of steps for the inner loop
1933      int oleftstep=0;  // the o variables refer to the outer loop
1934      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1935      int onumsteps=1;
1936      
1937      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1938      bool RES=(rightExp && rightScalar);
1939      bool LS=(!leftExp && leftScalar); // left is a single scalar
1940      bool RS=(!rightExp && rightScalar);
1941      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1942      bool RN=(!rightExp && !rightScalar);
1943      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1944      bool REN=(rightExp && !rightScalar);
1945    
1946      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1947      {
1948        chunksize=m_left->getNumDPPSample()*leftsize;
1949        leftstep=0;
1950        rightstep=0;
1951        numsteps=1;
1952      }
1953      else if (LES || RES)
1954      {
1955        chunksize=1;
1956        if (LES)        // left is an expanded scalar
1957        {
1958            if (RS)
1959            {
1960               leftstep=1;
1961               rightstep=0;
1962               numsteps=m_left->getNumDPPSample();
1963            }
1964            else        // RN or REN
1965            {
1966               leftstep=0;
1967               oleftstep=1;
1968               rightstep=1;
1969               orightstep=(RN ? -(int)rightsize : 0);
1970               numsteps=rightsize;
1971               onumsteps=m_left->getNumDPPSample();
1972            }
1973        }
1974        else        // right is an expanded scalar
1975        {
1976            if (LS)
1977            {
1978               rightstep=1;
1979               leftstep=0;
1980               numsteps=m_right->getNumDPPSample();
1981            }
1982            else
1983            {
1984               rightstep=0;
1985               orightstep=1;
1986               leftstep=1;
1987               oleftstep=(LN ? -(int)leftsize : 0);
1988               numsteps=leftsize;
1989               onumsteps=m_right->getNumDPPSample();
1990            }
1991        }
1992      }
1993      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1994      {
1995        if (LEN)    // and Right will be a single value
1996        {
1997            chunksize=rightsize;
1998            leftstep=rightsize;
1999            rightstep=0;
2000            numsteps=m_left->getNumDPPSample();
2001            if (RS)
2002            {
2003               numsteps*=leftsize;
2004            }
2005        }
2006        else    // REN
2007        {
2008            chunksize=leftsize;
2009            rightstep=leftsize;
2010            leftstep=0;
2011            numsteps=m_right->getNumDPPSample();
2012            if (LS)
2013            {
2014               numsteps*=rightsize;
2015            }
2016        }
2017      }
2018    
2019      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
2020        // Get the values of sub-expressions
2021      const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
2022      const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
2023    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
2024    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
2025    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
2026    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
2027    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
2028    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
2029    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
2030    
2031    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
2032    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
2033    
2034    
2035      roffset=m_samplesize*tid;
2036      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
2037      switch(m_op)
2038      {
2039        case ADD:
2040            PROC_OP(NO_ARG,plus<double>());
2041        break;
2042        case SUB:
2043        PROC_OP(NO_ARG,minus<double>());
2044        break;
2045        case MUL:
2046        PROC_OP(NO_ARG,multiplies<double>());
2047        break;
2048        case DIV:
2049        PROC_OP(NO_ARG,divides<double>());
2050        break;
2051        case POW:
2052           PROC_OP(double (double,double),::pow);
2053        break;
2054        default:
2055        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
2056      }
2057    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
2058      return &m_samples;
2059    }
2060    
2061    
2062    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
2063    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
2064    // unlike the other resolve helpers, we must treat these datapoints separately.
2065    const DataTypes::ValueType*
2066    DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
2067    {
2068    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
2069    
2070      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
2071        // first work out which of the children are expanded
2072      bool leftExp=(m_left->m_readytype=='E');
2073      bool rightExp=(m_right->m_readytype=='E');
2074      int steps=getNumDPPSample();
2075      int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
2076      int rightStep=(rightExp?m_right->getNoValues() : 0);
2077    
2078      int resultStep=getNoValues();
2079      roffset=m_samplesize*tid;
2080      size_t offset=roffset;
2081    
2082      const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
2083    
2084      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
2085    
2086    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
2087    cout << getNoValues() << endl;)
2088    
2089    
2090    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
2091    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
2092    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
2093    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
2094    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
2095    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
2096    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
2097    
2098      double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
2099      switch(m_op)
2100      {
2101        case PROD:
2102        for (int i=0;i<steps;++i,resultp+=resultStep)
2103        {
2104              const double *ptr_0 = &((*left)[lroffset]);
2105              const double *ptr_1 = &((*right)[rroffset]);
2106    
2107    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
2108    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
2109    
2110              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
2111    
2112          lroffset+=leftStep;
2113          rroffset+=rightStep;
2114        }
2115        break;
2116        default:
2117        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
2118      }
2119      roffset=offset;
2120      return &m_samples;
2121    }
2122    #endif //LAZY_NODE_STORAGE
2123    
2124  /*  /*
2125    \brief Compute the value of the expression for the given sample.    \brief Compute the value of the expression for the given sample.
# Line 1304  LAZYDEBUG(cout << "Resolve sample " << t Line 2148  LAZYDEBUG(cout << "Resolve sample " << t
2148    }    }
2149    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
2150    {    {
2151      const ValueType& vec=m_id->getVector();      const ValueType& vec=m_id->getVectorRO();
2152      if (m_readytype=='C')      if (m_readytype=='C')
2153      {      {
2154      roffset=0;      roffset=0;
# Line 1327  LAZYDEBUG(cout << "Finish  sample " << t Line 2171  LAZYDEBUG(cout << "Finish  sample " << t
2171    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
2172    case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);    case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
2173    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);    case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
2174      case G_NP1OUT_2P: return resolveNP1OUT_2P(v, offset, sampleNo, roffset);
2175    default:    default:
2176      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
2177    }    }
2178    
2179  }  }
2180    
2181    const DataTypes::ValueType*
2182    DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
2183    {
2184    #ifdef _OPENMP
2185        int tid=omp_get_thread_num();
2186    #else
2187        int tid=0;
2188    #endif
2189        return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
2190    }
2191    
2192    
2193    // This needs to do the work of the identity constructor
2194    void
2195    DataLazy::resolveToIdentity()
2196    {
2197       if (m_op==IDENTITY)
2198        return;
2199    #ifndef LAZY_NODE_STORAGE
2200       DataReady_ptr p=resolveVectorWorker();
2201    #else
2202       DataReady_ptr p=resolveNodeWorker();
2203    #endif
2204       makeIdentity(p);
2205    }
2206    
2207    void DataLazy::makeIdentity(const DataReady_ptr& p)
2208    {
2209       m_op=IDENTITY;
2210       m_axis_offset=0;
2211       m_transpose=0;
2212       m_SL=m_SM=m_SR=0;
2213       m_children=m_height=0;
2214       m_id=p;
2215       if(p->isConstant()) {m_readytype='C';}
2216       else if(p->isExpanded()) {m_readytype='E';}
2217       else if (p->isTagged()) {m_readytype='T';}
2218       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
2219       m_buffsRequired=1;
2220       m_samplesize=p->getNumDPPSample()*p->getNoValues();
2221       m_maxsamplesize=m_samplesize;
2222       m_left.reset();
2223       m_right.reset();
2224    }
2225    
2226    
2227    DataReady_ptr
2228    DataLazy::resolve()
2229    {
2230        resolveToIdentity();
2231        return m_id;
2232    }
2233    
2234    #ifdef LAZY_NODE_STORAGE
2235    
2236    // This version of resolve uses storage in each node to hold results
2237    DataReady_ptr
2238    DataLazy::resolveNodeWorker()
2239    {
2240      if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
2241      {
2242        collapse();
2243      }
2244      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
2245      {
2246        return m_id;
2247      }
2248        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
2249      DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
2250      ValueType& resvec=result->getVectorRW();
2251      DataReady_ptr resptr=DataReady_ptr(result);
2252    
2253      int sample;
2254      int totalsamples=getNumSamples();
2255      const ValueType* res=0;   // Storage for answer
2256    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
2257      #pragma omp parallel for private(sample,res) schedule(static)
2258      for (sample=0;sample<totalsamples;++sample)
2259      {
2260        size_t roffset=0;
2261    #ifdef _OPENMP
2262        res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
2263    #else
2264        res=resolveNodeSample(0,sample,roffset);
2265    #endif
2266    LAZYDEBUG(cout << "Sample #" << sample << endl;)
2267    LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
2268        DataVector::size_type outoffset=result->getPointOffset(sample,0);
2269        memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
2270      }
2271      return resptr;
2272    }
2273    
2274    #endif // LAZY_NODE_STORAGE
2275    
2276  // To simplify the memory management, all threads operate on one large vector, rather than one each.  // To simplify the memory management, all threads operate on one large vector, rather than one each.
2277  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
2278  DataReady_ptr  DataReady_ptr
2279  DataLazy::resolve()  DataLazy::resolveVectorWorker()
2280  {  {
2281    
2282  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
2283  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
   
2284    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
2285    {    {
2286      collapse();      collapse();
# Line 1356  LAZYDEBUG(cout << "Buffers=" << m_buffsR Line 2294  LAZYDEBUG(cout << "Buffers=" << m_buffsR
2294      // storage to evaluate its expression      // storage to evaluate its expression
2295    int numthreads=1;    int numthreads=1;
2296  #ifdef _OPENMP  #ifdef _OPENMP
2297    numthreads=getNumberOfThreads();    numthreads=omp_get_max_threads();
2298  #endif  #endif
2299    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
2300  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
2301    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
2302    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
2303    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
2304    int sample;    int sample;
2305    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
# Line 1372  LAZYDEBUG(cout << "Total number of sampl Line 2310  LAZYDEBUG(cout << "Total number of sampl
2310    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
2311    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
2312    {    {
       if (sample==0)  {ENABLEDEBUG}  
   
 //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}  
2313  LAZYDEBUG(cout << "################################# " << sample << endl;)  LAZYDEBUG(cout << "################################# " << sample << endl;)
2314  #ifdef _OPENMP  #ifdef _OPENMP
2315      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
# Line 1392  LAZYDEBUG(cerr << "outoffset=" << outoff Line 2327  LAZYDEBUG(cerr << "outoffset=" << outoff
2327      }      }
2328  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
2329  LAZYDEBUG(cerr << "*********************************" << endl;)  LAZYDEBUG(cerr << "*********************************" << endl;)
     DISABLEDEBUG  
2330    }    }
2331    return resptr;    return resptr;
2332  }  }
# Line 1410  DataLazy::toString() const Line 2344  DataLazy::toString() const
2344  void  void
2345  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
2346  {  {
2347    //    oss << "[" << m_children <<";"<<m_height <<"]";
2348    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
2349    {    {
2350    case G_IDENTITY:    case G_IDENTITY:
# Line 1453  DataLazy::intoString(ostringstream& oss) Line 2388  DataLazy::intoString(ostringstream& oss)
2388      m_right->intoString(oss);      m_right->intoString(oss);
2389      oss << ')';      oss << ')';
2390      break;      break;
2391      case G_NP1OUT_2P:
2392        oss << opToString(m_op) << '(';
2393        m_left->intoString(oss);
2394        oss << ", " << m_axis_offset << ", " << m_transpose;
2395        oss << ')';
2396        break;
2397    default:    default:
2398      oss << "UNKNOWN";      oss << "UNKNOWN";
2399    }    }
# Line 1549  DataLazy::getPointOffset(int sampleNo, Line 2490  DataLazy::getPointOffset(int sampleNo,
2490    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).");
2491  }  }
2492    
2493  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
2494  // 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.  
2495  void  void
2496  DataLazy::setToZero()  DataLazy::setToZero()
2497  {  {
2498    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
2499    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2500    m_op=IDENTITY;  //   m_op=IDENTITY;
2501    m_right.reset();    //   m_right.reset();  
2502    m_left.reset();  //   m_left.reset();
2503    m_readytype='C';  //   m_readytype='C';
2504    m_buffsRequired=1;  //   m_buffsRequired=1;
2505    
2506      privdebug=privdebug;  // to stop the compiler complaining about unused privdebug
2507      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2508    }
2509    
2510    bool
2511    DataLazy::actsExpanded() const
2512    {
2513        return (m_readytype=='E');
2514  }  }
2515    
2516  }   // end namespace  }   // end namespace

Legend:
Removed from v.2166  
changed lines
  Added in v.2514

  ViewVC Help
Powered by ViewVC 1.1.26