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

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

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

revision 2147 by jfenwick, Wed Dec 10 04:41:26 2008 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 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 16  Line 16 
16  #ifdef USE_NETCDF  #ifdef USE_NETCDF
17  #include <netcdfcpp.h>  #include <netcdfcpp.h>
18  #endif  #endif
19  #ifdef PASO_MPI  #ifdef ESYS_MPI
20  #include <mpi.h>  #include <mpi.h>
21  #endif  #endif
22  #ifdef _OPENMP  #ifdef _OPENMP
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31  // #define LAZYDEBUG(X) X;  #include "EscriptParams.h"
32    
33    #include <iomanip>      // for some fancy formatting in debug
34    
35    // #define LAZYDEBUG(X) if (privdebug){X;}
36  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
37    namespace
38    {
39    bool privdebug=false;
40    
41    #define ENABLEDEBUG privdebug=true;
42    #define DISABLEDEBUG privdebug=false;
43    }
44    
45    // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl;resolveToIdentity();}
46    
47    // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "SIZE LIMIT EXCEEDED " << m_height << endl;resolveToIdentity();}
48    
49    
50    #define SIZELIMIT if (m_height>escript::escriptParams.getTOO_MANY_LEVELS())  {if (escript::escriptParams.getLAZY_VERBOSE()){cerr << "SIZE LIMIT EXCEEDED height=" << m_height << endl;}resolveToIdentity();}
51    
52  /*  /*
53  How does DataLazy work?  How does DataLazy work?
# Line 93  namespace escript Line 111  namespace escript
111  namespace  namespace
112  {  {
113    
114    
115    // enabling this will print out when ever the maximum stacksize used by resolve increases
116    // it assumes _OPENMP is also in use
117    //#define LAZY_STACK_PROF
118    
119    
120    
121    #ifndef _OPENMP
122      #ifdef LAZY_STACK_PROF
123      #undef LAZY_STACK_PROF
124      #endif
125    #endif
126    
127    
128    #ifdef LAZY_STACK_PROF
129    std::vector<void*> stackstart(getNumberOfThreads());
130    std::vector<void*> stackend(getNumberOfThreads());
131    size_t maxstackuse=0;
132    #endif
133    
134  enum ES_opgroup  enum ES_opgroup
135  {  {
136     G_UNKNOWN,     G_UNKNOWN,
# Line 102  enum ES_opgroup Line 140  enum ES_opgroup
140     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
141     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
142     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
143     G_TENSORPROD     // general tensor product     G_TENSORPROD,    // general tensor product
144       G_NP1OUT_2P,     // non-pointwise op with one output requiring two params
145       G_REDUCTION,     // non-pointwise unary op with a scalar output
146       G_CONDEVAL
147  };  };
148    
149    
# Line 116  string ES_opstrings[]={"UNKNOWN","IDENTI Line 157  string ES_opstrings[]={"UNKNOWN","IDENTI
157              "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",
158              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
159              "prod",              "prod",
160              "transpose", "trace"};              "transpose", "trace",
161  int ES_opcount=40;              "swapaxes",
162                "minval", "maxval",
163                "condEval"};
164    int ES_opcount=44;
165  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,
166              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
167              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
# Line 126  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT Line 170  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT
170              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
171              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
172              G_TENSORPROD,              G_TENSORPROD,
173              G_NP1OUT_P, G_NP1OUT_P};              G_NP1OUT_P, G_NP1OUT_P,
174                G_NP1OUT_2P,
175                G_REDUCTION, G_REDUCTION,
176                G_CONDEVAL};
177  inline  inline
178  ES_opgroup  ES_opgroup
179  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 171  resultShape(DataAbstract_ptr left, DataA Line 218  resultShape(DataAbstract_ptr left, DataA
218        {        {
219          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.");
220        }        }
221    
222        if (left->getRank()==0)   // we need to allow scalar * anything        if (left->getRank()==0)   // we need to allow scalar * anything
223        {        {
224          return right->getShape();          return right->getShape();
# Line 187  resultShape(DataAbstract_ptr left, DataA Line 235  resultShape(DataAbstract_ptr left, DataA
235  // return the shape for "op left"  // return the shape for "op left"
236    
237  DataTypes::ShapeType  DataTypes::ShapeType
238  resultShape(DataAbstract_ptr left, ES_optype op)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
239  {  {
240      switch(op)      switch(op)
241      {      {
242          case TRANS:          case TRANS:
243          return left->getShape();         {            // for the scoping of variables
244            const DataTypes::ShapeType& s=left->getShape();
245            DataTypes::ShapeType sh;
246            int rank=left->getRank();
247            if (axis_offset<0 || axis_offset>rank)
248            {
249                stringstream e;
250                e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
251                throw DataException(e.str());
252            }
253            for (int i=0; i<rank; i++)
254            {
255               int index = (axis_offset+i)%rank;
256               sh.push_back(s[index]); // Append to new shape
257            }
258            return sh;
259           }
260      break;      break;
261      case TRACE:      case TRACE:
262          return DataTypes::scalarShape;         {
263            int rank=left->getRank();
264            if (rank<2)
265            {
266               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
267            }
268            if ((axis_offset>rank-2) || (axis_offset<0))
269            {
270               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
271            }
272            if (rank==2)
273            {
274               return DataTypes::scalarShape;
275            }
276            else if (rank==3)
277            {
278               DataTypes::ShapeType sh;
279                   if (axis_offset==0)
280               {
281                    sh.push_back(left->getShape()[2]);
282                   }
283                   else     // offset==1
284               {
285                sh.push_back(left->getShape()[0]);
286                   }
287               return sh;
288            }
289            else if (rank==4)
290            {
291               DataTypes::ShapeType sh;
292               const DataTypes::ShapeType& s=left->getShape();
293                   if (axis_offset==0)
294               {
295                    sh.push_back(s[2]);
296                    sh.push_back(s[3]);
297                   }
298                   else if (axis_offset==1)
299               {
300                    sh.push_back(s[0]);
301                    sh.push_back(s[3]);
302                   }
303               else     // offset==2
304               {
305                sh.push_back(s[0]);
306                sh.push_back(s[1]);
307               }
308               return sh;
309            }
310            else        // unknown rank
311            {
312               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
313            }
314           }
315      break;      break;
316          default:          default:
317      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)+".");
318      }      }
319  }  }
320    
321    DataTypes::ShapeType
322    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
323    {
324         // This code taken from the Data.cpp swapaxes() method
325         // Some of the checks are probably redundant here
326         int axis0_tmp,axis1_tmp;
327         const DataTypes::ShapeType& s=left->getShape();
328         DataTypes::ShapeType out_shape;
329         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
330         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
331         int rank=left->getRank();
332         if (rank<2) {
333            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
334         }
335         if (axis0<0 || axis0>rank-1) {
336            stringstream e;
337            e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
338            throw DataException(e.str());
339         }
340         if (axis1<0 || axis1>rank-1) {
341            stringstream e;
342            e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
343            throw DataException(e.str());
344         }
345         if (axis0 == axis1) {
346             throw DataException("Error - Data::swapaxes: axis indices must be different.");
347         }
348         if (axis0 > axis1) {
349             axis0_tmp=axis1;
350             axis1_tmp=axis0;
351         } else {
352             axis0_tmp=axis0;
353             axis1_tmp=axis1;
354         }
355         for (int i=0; i<rank; i++) {
356           if (i == axis0_tmp) {
357              out_shape.push_back(s[axis1_tmp]);
358           } else if (i == axis1_tmp) {
359              out_shape.push_back(s[axis0_tmp]);
360           } else {
361              out_shape.push_back(s[i]);
362           }
363         }
364        return out_shape;
365    }
366    
367    
368  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
369  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
370  // the majority of this code is copy pasted from C_General_Tensor_Product  // the majority of this code is copy pasted from C_General_Tensor_Product
# Line 266  GTPShape(DataAbstract_ptr left, DataAbst Line 429  GTPShape(DataAbstract_ptr left, DataAbst
429    return shape2;    return shape2;
430  }  }
431    
   
 // determine the number of points in the result of "left op right"  
 // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here  
 // size_t  
 // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
 // {  
 //    switch (getOpgroup(op))  
 //    {  
 //    case G_BINARY: return left->getLength();  
 //    case G_UNARY: return left->getLength();  
 //    case G_NP1OUT: return left->getLength();  
 //    default:  
 //  throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");  
 //    }  
 // }  
   
 // determine the number of samples requires to evaluate an expression combining left and right  
 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.  
 // The same goes for G_TENSORPROD  
 int  
 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  
 {  
    switch(getOpgroup(op))  
    {  
    case G_IDENTITY: return 1;  
    case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);  
    case G_UNARY:  
    case G_UNARY_P: return max(left->getBuffsRequired(),1);  
    case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);  
    case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);  
    case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);  
    default:  
     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");  
    }  
 }  
   
   
432  }   // end anonymous namespace  }   // end anonymous namespace
433    
434    
# Line 318  opToString(ES_optype op) Line 444  opToString(ES_optype op)
444    return ES_opstrings[op];    return ES_opstrings[op];
445  }  }
446    
447    void DataLazy::LazyNodeSetup()
448    {
449    #ifdef _OPENMP
450        int numthreads=omp_get_max_threads();
451        m_samples.resize(numthreads*m_samplesize);
452        m_sampleids=new int[numthreads];
453        for (int i=0;i<numthreads;++i)
454        {
455            m_sampleids[i]=-1;  
456        }
457    #else
458        m_samples.resize(m_samplesize);
459        m_sampleids=new int[1];
460        m_sampleids[0]=-1;
461    #endif  // _OPENMP
462    }
463    
464    
465    // Creates an identity node
466  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
467      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
468      m_op(IDENTITY),      ,m_sampleids(0),
469      m_axis_offset(0),      m_samples(1)
     m_transpose(0),  
     m_SL(0), m_SM(0), m_SR(0)  
470  {  {
471     if (p->isLazy())     if (p->isLazy())
472     {     {
# Line 335  DataLazy::DataLazy(DataAbstract_ptr p) Line 477  DataLazy::DataLazy(DataAbstract_ptr p)
477     }     }
478     else     else
479     {     {
480      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
481      if(p->isConstant()) {m_readytype='C';}      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
482      else if(p->isExpanded()) {m_readytype='E';}      makeIdentity(dr);
483      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.");}  
484     }     }
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_maxsamplesize=m_samplesize;  
485  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
486  }  }
487    
   
   
   
488  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
489      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
490      m_op(op),      m_op(op),
491      m_axis_offset(0),      m_axis_offset(0),
492      m_transpose(0),      m_transpose(0),
493      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
494  {  {
495     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
496     {     {
497      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.");
498     }     }
# Line 373  DataLazy::DataLazy(DataAbstract_ptr left Line 508  DataLazy::DataLazy(DataAbstract_ptr left
508     }     }
509     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
510     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
511     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
512     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
513       m_height=m_left->m_height+1;
514       LazyNodeSetup();
515       SIZELIMIT
516  }  }
517    
518    
# Line 385  DataLazy::DataLazy(DataAbstract_ptr left Line 522  DataLazy::DataLazy(DataAbstract_ptr left
522      m_op(op),      m_op(op),
523      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
524  {  {
525    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
526     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
527     {     {
528      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
# Line 401  DataLazy::DataLazy(DataAbstract_ptr left Line 539  DataLazy::DataLazy(DataAbstract_ptr left
539     {     {
540      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
541      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
542    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
543     }     }
544     left->operandCheck(*right);     left->operandCheck(*right);
545    
546     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
547     {     {
548      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
549    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
550     }     }
551     else     else
552     {     {
553      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
554    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
555     }     }
556     if (right->isLazy())     if (right->isLazy())
557     {     {
558      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
559    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
560     }     }
561     else     else
562     {     {
563      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
564    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
565     }     }
566     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
567     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 435  DataLazy::DataLazy(DataAbstract_ptr left Line 578  DataLazy::DataLazy(DataAbstract_ptr left
578      m_readytype='C';      m_readytype='C';
579     }     }
580     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
581     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_children=m_left->m_children+m_right->m_children+2;
582     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
583       LazyNodeSetup();
584       SIZELIMIT
585  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
586  }  }
587    
# Line 466  DataLazy::DataLazy(DataAbstract_ptr left Line 611  DataLazy::DataLazy(DataAbstract_ptr left
611      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
612      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
613     }     }
614     left->operandCheck(*right);  //    left->operandCheck(*right);
615    
616     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
617     {     {
# Line 499  DataLazy::DataLazy(DataAbstract_ptr left Line 644  DataLazy::DataLazy(DataAbstract_ptr left
644      m_readytype='C';      m_readytype='C';
645     }     }
646     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
647     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());       m_children=m_left->m_children+m_right->m_children+2;
648     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_height=max(m_left->m_height,m_right->m_height)+1;
649       LazyNodeSetup();
650       SIZELIMIT
651  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
652  }  }
653    
654    
655  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
656      : parent(left->getFunctionSpace(), resultShape(left,op)),      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
657      m_op(op),      m_op(op),
658      m_axis_offset(axis_offset),      m_axis_offset(axis_offset),
659      m_transpose(0),      m_transpose(0),
# Line 527  DataLazy::DataLazy(DataAbstract_ptr left Line 674  DataLazy::DataLazy(DataAbstract_ptr left
674     }     }
675     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
676     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
677     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
678     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
679       m_height=m_left->m_height+1;
680       LazyNodeSetup();
681       SIZELIMIT
682  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
683  }  }
684    
# Line 555  DataLazy::DataLazy(DataAbstract_ptr left Line 704  DataLazy::DataLazy(DataAbstract_ptr left
704     }     }
705     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
706     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
707     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
708     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_children=m_left->m_children+1;
709       m_height=m_left->m_height+1;
710       LazyNodeSetup();
711       SIZELIMIT
712  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
713  }  }
714    
715  DataLazy::~DataLazy()  
716    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
717        : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
718        m_op(op),
719        m_axis_offset(axis0),
720        m_transpose(axis1),
721        m_tol(0)
722  {  {
723       if ((getOpgroup(op)!=G_NP1OUT_2P))
724       {
725        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
726       }
727       DataLazy_ptr lleft;
728       if (!left->isLazy())
729       {
730        lleft=DataLazy_ptr(new DataLazy(left));
731       }
732       else
733       {
734        lleft=dynamic_pointer_cast<DataLazy>(left);
735       }
736       m_readytype=lleft->m_readytype;
737       m_left=lleft;
738       m_samplesize=getNumDPPSample()*getNoValues();
739       m_children=m_left->m_children+1;
740       m_height=m_left->m_height+1;
741       LazyNodeSetup();
742       SIZELIMIT
743    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
744  }  }
745    
746    
747  int  namespace
 DataLazy::getBuffsRequired() const  
748  {  {
749      return m_buffsRequired;  
750        inline int max3(int a, int b, int c)
751        {
752        int t=(a>b?a:b);
753        return (t>c?t:c);
754    
755        }
756  }  }
757    
758    DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
759        : parent(left->getFunctionSpace(), left->getShape()),
760        m_op(CONDEVAL),
761        m_axis_offset(0),
762        m_transpose(0),
763        m_tol(0)
764    {
765    
766  size_t     DataLazy_ptr lmask;
767  DataLazy::getMaxSampleSize() const     DataLazy_ptr lleft;
768       DataLazy_ptr lright;
769       if (!mask->isLazy())
770       {
771        lmask=DataLazy_ptr(new DataLazy(mask));
772       }
773       else
774       {
775        lmask=dynamic_pointer_cast<DataLazy>(mask);
776       }
777       if (!left->isLazy())
778       {
779        lleft=DataLazy_ptr(new DataLazy(left));
780       }
781       else
782       {
783        lleft=dynamic_pointer_cast<DataLazy>(left);
784       }
785       if (!right->isLazy())
786       {
787        lright=DataLazy_ptr(new DataLazy(right));
788       }
789       else
790       {
791        lright=dynamic_pointer_cast<DataLazy>(right);
792       }
793       m_readytype=lmask->m_readytype;
794       if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
795       {
796        throw DataException("Programmer Error - condEval arguments must have the same readytype");
797       }
798       m_left=lleft;
799       m_right=lright;
800       m_mask=lmask;
801       m_samplesize=getNumDPPSample()*getNoValues();
802       m_children=m_left->m_children+m_right->m_children+m_mask->m_children+1;
803       m_height=max3(m_left->m_height,m_right->m_height,m_mask->m_height)+1;
804       LazyNodeSetup();
805       SIZELIMIT
806    LAZYDEBUG(cout << "(8)Lazy created with " << m_samplesize << endl;)
807    }
808    
809    
810    
811    DataLazy::~DataLazy()
812  {  {
813      return m_maxsamplesize;     delete[] m_sampleids;
814  }  }
815    
816    
817  /*  /*
818    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
819    This does the work for the collapse method.    This does the work for the collapse method.
# Line 718  DataLazy::collapseToReady() Line 953  DataLazy::collapseToReady()
953      case TRACE:      case TRACE:
954      result=left.trace(m_axis_offset);      result=left.trace(m_axis_offset);
955      break;      break;
956        case SWAP:
957        result=left.swapaxes(m_axis_offset, m_transpose);
958        break;
959        case MINVAL:
960        result=left.minval();
961        break;
962        case MAXVAL:
963        result=left.minval();
964        break;
965      default:      default:
966      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)+".");
967    }    }
# Line 745  DataLazy::collapse() Line 989  DataLazy::collapse()
989    m_op=IDENTITY;    m_op=IDENTITY;
990  }  }
991    
992  /*  
993    \brief Compute the value of the expression (unary operation) for the given sample.  
994    \return Vector which stores the value of the subexpression for the given sample.  
995    \param v A vector to store intermediate results.  
996    \param offset Index in v to begin storing results.  
997    \param sampleNo Sample number to evaluate.  #define PROC_OP(TYPE,X)                               \
998    \param roffset (output parameter) the offset in the return vector where the result begins.      for (int j=0;j<onumsteps;++j)\
999        {\
1000    The return value will be an existing vector so do not deallocate it.        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1001    If the result is stored in v it should be stored at the offset given.        { \
1002    Everything from offset to the end of v should be considered available for this method to use.  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1003  */  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1004  DataTypes::ValueType*           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1005  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1006             lroffset+=leftstep; \
1007             rroffset+=rightstep; \
1008          }\
1009          lroffset+=oleftstep;\
1010          rroffset+=orightstep;\
1011        }
1012    
1013    
1014    // The result will be stored in m_samples
1015    // The return value is a pointer to the DataVector, offset is the offset within the return value
1016    const DataTypes::ValueType*
1017    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
1018    {
1019    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1020        // collapse so we have a 'E' node or an IDENTITY for some other type
1021      if (m_readytype!='E' && m_op!=IDENTITY)
1022      {
1023        collapse();
1024      }
1025      if (m_op==IDENTITY)  
1026      {
1027        const ValueType& vec=m_id->getVectorRO();
1028        roffset=m_id->getPointOffset(sampleNo, 0);
1029    #ifdef LAZY_STACK_PROF
1030    int x;
1031    if (&x<stackend[omp_get_thread_num()])
1032    {
1033           stackend[omp_get_thread_num()]=&x;
1034    }
1035    #endif
1036        return &(vec);
1037      }
1038      if (m_readytype!='E')
1039      {
1040        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1041      }
1042      if (m_sampleids[tid]==sampleNo)
1043      {
1044        roffset=tid*m_samplesize;
1045        return &(m_samples);        // sample is already resolved
1046      }
1047      m_sampleids[tid]=sampleNo;
1048    
1049      switch (getOpgroup(m_op))
1050      {
1051      case G_UNARY:
1052      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1053      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1054      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1055      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1056      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1057      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1058      case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
1059      case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
1060      default:
1061        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1062      }
1063    }
1064    
1065    const DataTypes::ValueType*
1066    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
1067  {  {
1068      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1069      // 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
1070      // processing single points.      // processing single points.
1071        // we will also know we won't get identity nodes
1072    if (m_readytype!='E')    if (m_readytype!='E')
1073    {    {
1074      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1075    }    }
1076    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);    if (m_op==IDENTITY)
1077    const double* left=&((*vleft)[roffset]);    {
1078    double* result=&(v[offset]);      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1079    roffset=offset;    }
1080      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1081      const double* left=&((*leftres)[roffset]);
1082      roffset=m_samplesize*tid;
1083      double* result=&(m_samples[roffset]);
1084    switch (m_op)    switch (m_op)
1085    {    {
1086      case SIN:        case SIN:  
# Line 880  DataLazy::resolveUnary(ValueType& v, siz Line 1190  DataLazy::resolveUnary(ValueType& v, siz
1190      default:      default:
1191      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1192    }    }
1193    return &v;    return &(m_samples);
1194  }  }
1195    
1196    
1197    const DataTypes::ValueType*
1198    DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)
1199    {
1200        // we assume that any collapsing has been done before we get here
1201        // since we only have one argument we don't need to think about only
1202        // processing single points.
1203        // we will also know we won't get identity nodes
1204      if (m_readytype!='E')
1205      {
1206        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1207      }
1208      if (m_op==IDENTITY)
1209      {
1210        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1211      }
1212      size_t loffset=0;
1213      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1214    
1215      roffset=m_samplesize*tid;
1216      unsigned int ndpps=getNumDPPSample();
1217      unsigned int psize=DataTypes::noValues(getShape());
1218      double* result=&(m_samples[roffset]);
1219      switch (m_op)
1220      {
1221        case MINVAL:
1222        {
1223          for (unsigned int z=0;z<ndpps;++z)
1224          {
1225            FMin op;
1226            *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1227            loffset+=psize;
1228            result++;
1229          }
1230        }
1231        break;
1232        case MAXVAL:
1233        {
1234          for (unsigned int z=0;z<ndpps;++z)
1235          {
1236          FMax op;
1237          *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1238          loffset+=psize;
1239          result++;
1240          }
1241        }
1242        break;
1243        default:
1244        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1245      }
1246      return &(m_samples);
1247    }
1248    
1249    const DataTypes::ValueType*
1250    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
 /*  
   \brief Compute the value of the expression (unary operation) for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
1251  {  {
1252      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1253      // 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
1254      // processing single points.      // processing single points.
1255    if (m_readytype!='E')    if (m_readytype!='E')
1256    {    {
1257      throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1258    }    }
1259      // since we can't write the result over the input, we need a result offset further along    if (m_op==IDENTITY)
1260    size_t subroffset=roffset+m_samplesize;    {
1261    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1262    roffset=offset;    }
1263      size_t subroffset;
1264      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1265      roffset=m_samplesize*tid;
1266      size_t loop=0;
1267      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1268      size_t step=getNoValues();
1269      size_t offset=roffset;
1270    switch (m_op)    switch (m_op)
1271    {    {
1272      case SYM:      case SYM:
1273      DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1274        {
1275            DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1276            subroffset+=step;
1277            offset+=step;
1278        }
1279      break;      break;
1280      case NSYM:      case NSYM:
1281      DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);      for (loop=0;loop<numsteps;++loop)
1282        {
1283            DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1284            subroffset+=step;
1285            offset+=step;
1286        }
1287      break;      break;
1288      default:      default:
1289      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1290    }    }
1291    return &v;    return &m_samples;
1292  }  }
1293    
1294  /*  const DataTypes::ValueType*
1295    \brief Compute the value of the expression (unary operation) for the given sample.  DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   The return value will be an existing vector so do not deallocate it.  
   If the result is stored in v it should be stored at the offset given.  
   Everything from offset to the end of v should be considered available for this method to use.  
 */  
 DataTypes::ValueType*  
 DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  
1296  {  {
1297      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1298      // 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
1299      // processing single points.      // processing single points.
1300    if (m_readytype!='E')    if (m_readytype!='E')
1301    {    {
1302      throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1303    }    }
1304      // since we can't write the result over the input, we need a result offset further along    if (m_op==IDENTITY)
1305    size_t subroffset=roffset+m_samplesize;    {
1306    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);      throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1307    roffset=offset;    }
1308      size_t subroffset;
1309      size_t offset;
1310      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1311      roffset=m_samplesize*tid;
1312      offset=roffset;
1313      size_t loop=0;
1314      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1315      size_t outstep=getNoValues();
1316      size_t instep=m_left->getNoValues();
1317    switch (m_op)    switch (m_op)
1318    {    {
1319      case TRACE:      case TRACE:
1320           DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
1321        {
1322                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1323            subroffset+=instep;
1324            offset+=outstep;
1325        }
1326      break;      break;
1327      case TRANS:      case TRANS:
1328           DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);      for (loop=0;loop<numsteps;++loop)
1329        {
1330                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1331            subroffset+=instep;
1332            offset+=outstep;
1333        }
1334      break;      break;
1335      default:      default:
1336      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1337    }    }
1338    return &v;    return &m_samples;
1339  }  }
1340    
1341    
1342  #define PROC_OP(TYPE,X)                               \  const DataTypes::ValueType*
1343      for (int i=0;i<steps;++i,resultp+=resultStep) \  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1344      { \  {
1345         tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    if (m_readytype!='E')
1346         lroffset+=leftStep; \    {
1347         rroffset+=rightStep; \      throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1348      }
1349      if (m_op==IDENTITY)
1350      {
1351        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1352      }
1353      size_t subroffset;
1354      size_t offset;
1355      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1356      roffset=m_samplesize*tid;
1357      offset=roffset;
1358      size_t loop=0;
1359      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1360      size_t outstep=getNoValues();
1361      size_t instep=m_left->getNoValues();
1362      switch (m_op)
1363      {
1364        case SWAP:
1365        for (loop=0;loop<numsteps;++loop)
1366        {
1367                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1368            subroffset+=instep;
1369            offset+=outstep;
1370      }      }
1371        break;
1372        default:
1373        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1374      }
1375      return &m_samples;
1376    }
1377    
1378    const DataTypes::ValueType*
1379    DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset)
1380    {
1381      if (m_readytype!='E')
1382      {
1383        throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1384      }
1385      if (m_op!=CONDEVAL)
1386      {
1387        throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1388      }
1389      size_t subroffset;
1390      size_t offset;
1391    
1392    
1393      const ValueType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1394      const ValueType* srcres=0;
1395      if ((*maskres)[subroffset]>0)
1396      {
1397        srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1398      }
1399      else
1400      {
1401        srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1402      }
1403    
1404      // Now we need to copy the result
1405    
1406      roffset=m_samplesize*tid;
1407      offset=roffset;
1408      for (int i=0;i<m_samplesize;++i)
1409      {
1410        m_samples[roffset+i]=(*srcres)[subroffset+i];  
1411      }
1412    
1413      return &m_samples;
1414    }
1415    
 /*  
   \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.  
 */  
1416  // 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
1417  // 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.
1418  // If both children are expanded, then we can process them in a single operation (we treat  // If both children are expanded, then we can process them in a single operation (we treat
# Line 998  DataLazy::resolveNP1OUT_P(ValueType& v, Line 1422  DataLazy::resolveNP1OUT_P(ValueType& v,
1422  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1423  // For example, 2+Vector.  // For example, 2+Vector.
1424  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1425  DataTypes::ValueType*  const DataTypes::ValueType*
1426  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1427  {  {
1428  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1429    
# Line 1013  LAZYDEBUG(cout << "Resolve binary: " << Line 1437  LAZYDEBUG(cout << "Resolve binary: " <<
1437    }    }
1438    bool leftScalar=(m_left->getRank()==0);    bool leftScalar=(m_left->getRank()==0);
1439    bool rightScalar=(m_right->getRank()==0);    bool rightScalar=(m_right->getRank()==0);
1440    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1441    int steps=(bigloops?1:getNumDPPSample());    {
1442    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint      throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1443    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops    }
1444    {    size_t leftsize=m_left->getNoValues();
1445      if (!leftScalar && !rightScalar)    size_t rightsize=m_right->getNoValues();
1446      {    size_t chunksize=1;           // how many doubles will be processed in one go
1447         throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");    int leftstep=0;       // how far should the left offset advance after each step
1448      }    int rightstep=0;
1449      steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());    int numsteps=0;       // total number of steps for the inner loop
1450      chunksize=1;    // for scalar    int oleftstep=0;  // the o variables refer to the outer loop
1451    }        int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1452    int leftStep=((leftExp && (!rightExp || rightScalar))? m_right->getNoValues() : 0);    int onumsteps=1;
1453    int rightStep=((rightExp && (!leftExp || leftScalar))? m_left->getNoValues() : 0);    
1454    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0    bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1455      bool RES=(rightExp && rightScalar);
1456      bool LS=(!leftExp && leftScalar); // left is a single scalar
1457      bool RS=(!rightExp && rightScalar);
1458      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1459      bool RN=(!rightExp && !rightScalar);
1460      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1461      bool REN=(rightExp && !rightScalar);
1462    
1463      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1464      {
1465        chunksize=m_left->getNumDPPSample()*leftsize;
1466        leftstep=0;
1467        rightstep=0;
1468        numsteps=1;
1469      }
1470      else if (LES || RES)
1471      {
1472        chunksize=1;
1473        if (LES)        // left is an expanded scalar
1474        {
1475            if (RS)
1476            {
1477               leftstep=1;
1478               rightstep=0;
1479               numsteps=m_left->getNumDPPSample();
1480            }
1481            else        // RN or REN
1482            {
1483               leftstep=0;
1484               oleftstep=1;
1485               rightstep=1;
1486               orightstep=(RN ? -(int)rightsize : 0);
1487               numsteps=rightsize;
1488               onumsteps=m_left->getNumDPPSample();
1489            }
1490        }
1491        else        // right is an expanded scalar
1492        {
1493            if (LS)
1494            {
1495               rightstep=1;
1496               leftstep=0;
1497               numsteps=m_right->getNumDPPSample();
1498            }
1499            else
1500            {
1501               rightstep=0;
1502               orightstep=1;
1503               leftstep=1;
1504               oleftstep=(LN ? -(int)leftsize : 0);
1505               numsteps=leftsize;
1506               onumsteps=m_right->getNumDPPSample();
1507            }
1508        }
1509      }
1510      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1511      {
1512        if (LEN)    // and Right will be a single value
1513        {
1514            chunksize=rightsize;
1515            leftstep=rightsize;
1516            rightstep=0;
1517            numsteps=m_left->getNumDPPSample();
1518            if (RS)
1519            {
1520               numsteps*=leftsize;
1521            }
1522        }
1523        else    // REN
1524        {
1525            chunksize=leftsize;
1526            rightstep=leftsize;
1527            leftstep=0;
1528            numsteps=m_right->getNumDPPSample();
1529            if (LS)
1530            {
1531               numsteps*=rightsize;
1532            }
1533        }
1534      }
1535    
1536      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1537      // Get the values of sub-expressions      // Get the values of sub-expressions
1538    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
1539    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1540      // the right child starts further along.  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1541    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;)
1542    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1543    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1544    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1545    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1546    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1547    
1548    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1549    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1550    
1551    
1552      roffset=m_samplesize*tid;
1553      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
1554    switch(m_op)    switch(m_op)
1555    {    {
1556      case ADD:      case ADD:
# Line 1053  LAZYDEBUG(cout << "Resolve binary: " << Line 1571  LAZYDEBUG(cout << "Resolve binary: " <<
1571      default:      default:
1572      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1573    }    }
1574    roffset=offset;    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1575    return &v;    return &m_samples;
1576  }  }
1577    
1578    
 /*  
   \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.  
 */  
1579  // 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
1580  // 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.
1581  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1582  DataTypes::ValueType*  const DataTypes::ValueType*
1583  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
1584  {  {
1585  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1586    
# Line 1083  LAZYDEBUG(cout << "Resolve TensorProduct Line 1589  LAZYDEBUG(cout << "Resolve TensorProduct
1589    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1590    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1591    int steps=getNumDPPSample();    int steps=getNumDPPSample();
1592    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1593    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1594    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0  
1595      // Get the values of sub-expressions (leave a gap of one sample for the result).    int resultStep=getNoValues();
1596    const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);    roffset=m_samplesize*tid;
1597    const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);    size_t offset=roffset;
1598    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved  
1599      const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1600    
1601      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1602    
1603    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1604    cout << getNoValues() << endl;)
1605    
1606    
1607    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1608    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1609    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1610    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1611    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1612    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1613    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1614    
1615      double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
1616    switch(m_op)    switch(m_op)
1617    {    {
1618      case PROD:      case PROD:
# Line 1097  LAZYDEBUG(cout << "Resolve TensorProduct Line 1620  LAZYDEBUG(cout << "Resolve TensorProduct
1620      {      {
1621            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1622            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1623    
1624    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1625    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1626    
1627            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);
1628    
1629        lroffset+=leftStep;        lroffset+=leftStep;
1630        rroffset+=rightStep;        rroffset+=rightStep;
1631      }      }
# Line 1106  LAZYDEBUG(cout << "Resolve TensorProduct Line 1634  LAZYDEBUG(cout << "Resolve TensorProduct
1634      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1635    }    }
1636    roffset=offset;    roffset=offset;
1637    return &v;    return &m_samples;
1638  }  }
1639    
1640    
1641    const DataTypes::ValueType*
1642    DataLazy::resolveSample(int sampleNo, size_t& roffset)
1643    {
1644    #ifdef _OPENMP
1645        int tid=omp_get_thread_num();
1646    #else
1647        int tid=0;
1648    #endif
1649    
1650  /*  #ifdef LAZY_STACK_PROF
1651    \brief Compute the value of the expression for the given sample.      stackstart[tid]=&tid;
1652    \return Vector which stores the value of the subexpression for the given sample.      stackend[tid]=&tid;
1653    \param v A vector to store intermediate results.      const DataTypes::ValueType* r=resolveNodeSample(tid, sampleNo, roffset);
1654    \param offset Index in v to begin storing results.      size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1655    \param sampleNo Sample number to evaluate.      #pragma omp critical
1656    \param roffset (output parameter) the offset in the return vector where the result begins.      if (d>maxstackuse)
1657        {
1658    cout << "Max resolve Stack use " << d << endl;
1659            maxstackuse=d;
1660        }
1661        return r;
1662    #else
1663        return resolveNodeSample(tid, sampleNo, roffset);
1664    #endif
1665    }
1666    
   The return value will be an existing vector so do not deallocate it.  
 */  
 // the vector and the offset are a place where the method could write its data if it wishes  
 // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
 // Hence the return value to indicate where the data is actually stored.  
 // Regardless, the storage should be assumed to be used, even if it isn't.  
1667    
1668  // the roffset is the offset within the returned vector where the data begins  // This needs to do the work of the identity constructor
1669  const DataTypes::ValueType*  void
1670  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveToIdentity()
1671  {  {
1672  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)     if (m_op==IDENTITY)
1673      // collapse so we have a 'E' node or an IDENTITY for some other type      return;
1674    if (m_readytype!='E' && m_op!=IDENTITY)     DataReady_ptr p=resolveNodeWorker();
1675       makeIdentity(p);
1676    }
1677    
1678    void DataLazy::makeIdentity(const DataReady_ptr& p)
1679    {
1680       m_op=IDENTITY;
1681       m_axis_offset=0;
1682       m_transpose=0;
1683       m_SL=m_SM=m_SR=0;
1684       m_children=m_height=0;
1685       m_id=p;
1686       if(p->isConstant()) {m_readytype='C';}
1687       else if(p->isExpanded()) {m_readytype='E';}
1688       else if (p->isTagged()) {m_readytype='T';}
1689       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1690       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1691       m_left.reset();
1692       m_right.reset();
1693    }
1694    
1695    
1696    DataReady_ptr
1697    DataLazy::resolve()
1698    {
1699        resolveToIdentity();
1700        return m_id;
1701    }
1702    
1703    
1704    /* This is really a static method but I think that caused problems in windows */
1705    void
1706    DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
1707    {
1708      if (dats.empty())
1709    {    {
1710      collapse();      return;
1711    }    }
1712    if (m_op==IDENTITY)      vector<DataLazy*> work;
1713      FunctionSpace fs=dats[0]->getFunctionSpace();
1714      bool match=true;
1715      for (int i=dats.size()-1;i>=0;--i)
1716    {    {
1717      const ValueType& vec=m_id->getVector();      if (dats[i]->m_readytype!='E')
1718      if (m_readytype=='C')      {
1719      {          dats[i]->collapse();
1720      roffset=0;      }
1721      return &(vec);      if (dats[i]->m_op!=IDENTITY)
1722      }      {
1723      roffset=m_id->getPointOffset(sampleNo, 0);          work.push_back(dats[i]);
1724      return &(vec);          if (fs!=dats[i]->getFunctionSpace())
1725            {
1726                match=false;
1727            }
1728        }
1729    }    }
1730    if (m_readytype!='E')    if (work.empty())
1731    {    {
1732      throw DataException("Programmer Error - Collapse did not produce an expanded node.");      return;     // no work to do
1733    }    }
1734    switch (getOpgroup(m_op))    if (match)    // all functionspaces match.  Yes I realise this is overly strict
1735      {     // it is possible that dats[0] is one of the objects which we discarded and
1736            // all the other functionspaces match.
1737        vector<DataExpanded*> dep;
1738        vector<ValueType*> vecs;
1739        for (int i=0;i<work.size();++i)
1740        {
1741            dep.push_back(new DataExpanded(fs,work[i]->getShape(), ValueType(work[i]->getNoValues())));
1742            vecs.push_back(&(dep[i]->getVectorRW()));
1743        }
1744        int totalsamples=work[0]->getNumSamples();
1745        const ValueType* res=0; // Storage for answer
1746        int sample;
1747        #pragma omp parallel private(sample, res)
1748        {
1749            size_t roffset=0;
1750            #pragma omp for schedule(static)
1751            for (sample=0;sample<totalsamples;++sample)
1752            {
1753            roffset=0;
1754            int j;
1755            for (j=work.size()-1;j>=0;--j)
1756            {
1757    #ifdef _OPENMP
1758                    res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1759    #else
1760                    res=work[j]->resolveNodeSample(0,sample,roffset);
1761    #endif
1762                    DataVector::size_type outoffset=dep[j]->getPointOffset(sample,0);
1763                    memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(DataVector::ElementType));
1764            }
1765            }
1766        }
1767        // Now we need to load the new results as identity ops into the lazy nodes
1768        for (int i=work.size()-1;i>=0;--i)
1769        {
1770            work[i]->makeIdentity(boost::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1771        }
1772      }
1773      else  // functionspaces do not match
1774    {    {
1775    case G_UNARY:      for (int i=0;i<work.size();++i)
1776    case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);      {
1777    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);          work[i]->resolveToIdentity();
1778    case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);      }
   case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);  
   case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);  
   default:  
     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
1779    }    }
1780  }  }
1781    
1782    
1783  // To simplify the memory management, all threads operate on one large vector, rather than one each.  
1784  // Each sample is evaluated independently and copied into the result DataExpanded.  // This version of resolve uses storage in each node to hold results
1785  DataReady_ptr  DataReady_ptr
1786  DataLazy::resolve()  DataLazy::resolveNodeWorker()
1787  {  {
   
 LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)  
 LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)  
   
1788    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
1789    {    {
1790      collapse();      collapse();
# Line 1183  LAZYDEBUG(cout << "Buffers=" << m_buffsR Line 1794  LAZYDEBUG(cout << "Buffers=" << m_buffsR
1794      return m_id;      return m_id;
1795    }    }
1796      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
   size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough  
     // storage to evaluate its expression  
   int numthreads=1;  
 #ifdef _OPENMP  
   numthreads=getNumberOfThreads();  
 #endif  
   ValueType v(numthreads*threadbuffersize);  
 LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)  
1797    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1798    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVectorRW();
1799    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1800    
1801    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1802    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1803    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Storage for answer
1804    size_t resoffset=0;       // where in the vector to find the answer  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1805    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel private(sample,res)
1806    for (sample=0;sample<totalsamples;++sample)    {
1807    {      size_t roffset=0;
1808  LAZYDEBUG(cout << "################################# " << sample << endl;)  #ifdef LAZY_STACK_PROF
1809        stackstart[omp_get_thread_num()]=&roffset;
1810        stackend[omp_get_thread_num()]=&roffset;
1811    #endif
1812        #pragma omp for schedule(static)
1813        for (sample=0;sample<totalsamples;++sample)
1814        {
1815            roffset=0;
1816  #ifdef _OPENMP  #ifdef _OPENMP
1817      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);              res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1818  #else  #else
1819      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.              res=resolveNodeSample(0,sample,roffset);
1820  #endif  #endif
1821  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1822      outoffset=result->getPointOffset(sample,0);  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1823  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)              DataVector::size_type outoffset=result->getPointOffset(sample,0);
1824      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));
1825      {      }
1826      resvec[outoffset]=(*res)[resoffset];    }
1827      }  #ifdef LAZY_STACK_PROF
1828  LAZYDEBUG(cerr << "*********************************" << endl;)    for (int i=0;i<getNumberOfThreads();++i)
1829      {
1830        size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1831    //  cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1832        if (r>maxstackuse)
1833        {
1834            maxstackuse=r;
1835        }
1836    }    }
1837      cout << "Max resolve Stack use=" << maxstackuse << endl;
1838    #endif
1839    return resptr;    return resptr;
1840  }  }
1841    
# Line 1224  std::string Line 1843  std::string
1843  DataLazy::toString() const  DataLazy::toString() const
1844  {  {
1845    ostringstream oss;    ostringstream oss;
1846    oss << "Lazy Data:";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1847    intoString(oss);    switch (escriptParams.getLAZY_STR_FMT())
1848      {
1849      case 1:   // tree format
1850        oss << endl;
1851        intoTreeString(oss,"");
1852        break;
1853      case 2:   // just the depth
1854        break;
1855      default:
1856        intoString(oss);
1857        break;
1858      }
1859    return oss.str();    return oss.str();
1860  }  }
1861    
# Line 1233  DataLazy::toString() const Line 1863  DataLazy::toString() const
1863  void  void
1864  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1865  {  {
1866    //    oss << "[" << m_children <<";"<<m_height <<"]";
1867    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1868    {    {
1869    case G_IDENTITY:    case G_IDENTITY:
# Line 1265  DataLazy::intoString(ostringstream& oss) Line 1896  DataLazy::intoString(ostringstream& oss)
1896    case G_UNARY_P:    case G_UNARY_P:
1897    case G_NP1OUT:    case G_NP1OUT:
1898    case G_NP1OUT_P:    case G_NP1OUT_P:
1899      case G_REDUCTION:
1900      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1901      m_left->intoString(oss);      m_left->intoString(oss);
1902      oss << ')';      oss << ')';
# Line 1276  DataLazy::intoString(ostringstream& oss) Line 1908  DataLazy::intoString(ostringstream& oss)
1908      m_right->intoString(oss);      m_right->intoString(oss);
1909      oss << ')';      oss << ')';
1910      break;      break;
1911      case G_NP1OUT_2P:
1912        oss << opToString(m_op) << '(';
1913        m_left->intoString(oss);
1914        oss << ", " << m_axis_offset << ", " << m_transpose;
1915        oss << ')';
1916        break;
1917      case G_CONDEVAL:
1918        oss << opToString(m_op)<< '(' ;
1919        m_mask->intoString(oss);
1920        oss << " ? ";
1921        m_left->intoString(oss);
1922        oss << " : ";
1923        m_right->intoString(oss);
1924        oss << ')';
1925        break;
1926    default:    default:
1927      oss << "UNKNOWN";      oss << "UNKNOWN";
1928    }    }
1929  }  }
1930    
1931    
1932    void
1933    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1934    {
1935      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1936      switch (getOpgroup(m_op))
1937      {
1938      case G_IDENTITY:
1939        if (m_id->isExpanded())
1940        {
1941           oss << "E";
1942        }
1943        else if (m_id->isTagged())
1944        {
1945          oss << "T";
1946        }
1947        else if (m_id->isConstant())
1948        {
1949          oss << "C";
1950        }
1951        else
1952        {
1953          oss << "?";
1954        }
1955        oss << '@' << m_id.get() << endl;
1956        break;
1957      case G_BINARY:
1958        oss << opToString(m_op) << endl;
1959        indent+='.';
1960        m_left->intoTreeString(oss, indent);
1961        m_right->intoTreeString(oss, indent);
1962        break;
1963      case G_UNARY:
1964      case G_UNARY_P:
1965      case G_NP1OUT:
1966      case G_NP1OUT_P:
1967      case G_REDUCTION:
1968        oss << opToString(m_op) << endl;
1969        indent+='.';
1970        m_left->intoTreeString(oss, indent);
1971        break;
1972      case G_TENSORPROD:
1973        oss << opToString(m_op) << endl;
1974        indent+='.';
1975        m_left->intoTreeString(oss, indent);
1976        m_right->intoTreeString(oss, indent);
1977        break;
1978      case G_NP1OUT_2P:
1979        oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1980        indent+='.';
1981        m_left->intoTreeString(oss, indent);
1982        break;
1983      default:
1984        oss << "UNKNOWN";
1985      }
1986    }
1987    
1988    
1989  DataAbstract*  DataAbstract*
1990  DataLazy::deepCopy()  DataLazy::deepCopy()
1991  {  {
1992    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1993    {    {
1994    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1995    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY:
1996      case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1997      case G_UNARY_P:   return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1998    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);
1999    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);
2000    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);
2001      case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
2002      case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2003    default:    default:
2004      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)+".");
2005    }    }
2006  }  }
2007    
2008    
2009    
2010  // There is no single, natural interpretation of getLength on DataLazy.  // There is no single, natural interpretation of getLength on DataLazy.
2011  // Instances of DataReady can look at the size of their vectors.  // Instances of DataReady can look at the size of their vectors.
2012  // For lazy though, it could be the size the data would be if it were resolved;  // For lazy though, it could be the size the data would be if it were resolved;
# Line 1372  DataLazy::getPointOffset(int sampleNo, Line 2082  DataLazy::getPointOffset(int sampleNo,
2082    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).");
2083  }  }
2084    
2085  // It would seem that DataTagged will need to be treated differently since even after setting all tags  
2086  // 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.  
2087  void  void
2088  DataLazy::setToZero()  DataLazy::setToZero()
2089  {  {
2090    DataTypes::ValueType v(getNoValues(),0);  //   DataTypes::ValueType v(getNoValues(),0);
2091    m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));  //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2092    m_op=IDENTITY;  //   m_op=IDENTITY;
2093    m_right.reset();    //   m_right.reset();  
2094    m_left.reset();  //   m_left.reset();
2095    m_readytype='C';  //   m_readytype='C';
2096    m_buffsRequired=1;  //   m_buffsRequired=1;
2097    
2098      privdebug=privdebug;  // to stop the compiler complaining about unused privdebug
2099      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2100    }
2101    
2102    bool
2103    DataLazy::actsExpanded() const
2104    {
2105        return (m_readytype=='E');
2106  }  }
2107    
2108  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26