/[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 2472 by jfenwick, Thu Jun 11 23:33:47 2009 UTC revision 4264 by jfenwick, Thu Feb 28 11:32:57 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    
17  #include "DataLazy.h"  #include "DataLazy.h"
18  #ifdef USE_NETCDF  #include "esysUtils/Esys_MPI.h"
 #include <netcdfcpp.h>  
 #endif  
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
19  #ifdef _OPENMP  #ifdef _OPENMP
20  #include <omp.h>  #include <omp.h>
21  #endif  #endif
# Line 30  Line 27 
27    
28  #include "EscriptParams.h"  #include "EscriptParams.h"
29    
30    #ifdef USE_NETCDF
31    #include <netcdfcpp.h>
32    #endif
33    
34  #include <iomanip>      // for some fancy formatting in debug  #include <iomanip>      // for some fancy formatting in debug
35    
36  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
# Line 42  bool privdebug=false; Line 43  bool privdebug=false;
43  #define DISABLEDEBUG privdebug=false;  #define DISABLEDEBUG privdebug=false;
44  }  }
45    
46  // #define SIZELIMIT  // #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();}
47  // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}  
48  //#define SIZELIMIT if ((m_height>7) || (m_children>127)) {resolveToIdentity();}  // #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();}
49    
50    
51  #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}  #define SIZELIMIT if (m_height>escript::escriptParams.getTOO_MANY_LEVELS())  {if (escript::escriptParams.getLAZY_VERBOSE()){cerr << "SIZE LIMIT EXCEEDED height=" << m_height << endl;}resolveToIdentity();}
52    
53  /*  /*
54  How does DataLazy work?  How does DataLazy work?
# Line 110  namespace escript Line 112  namespace escript
112  namespace  namespace
113  {  {
114    
115    
116    // enabling this will print out when ever the maximum stacksize used by resolve increases
117    // it assumes _OPENMP is also in use
118    //#define LAZY_STACK_PROF
119    
120    
121    
122    #ifndef _OPENMP
123      #ifdef LAZY_STACK_PROF
124      #undef LAZY_STACK_PROF
125      #endif
126    #endif
127    
128    
129    #ifdef LAZY_STACK_PROF
130    std::vector<void*> stackstart(getNumberOfThreads());
131    std::vector<void*> stackend(getNumberOfThreads());
132    size_t maxstackuse=0;
133    #endif
134    
135  enum ES_opgroup  enum ES_opgroup
136  {  {
137     G_UNKNOWN,     G_UNKNOWN,
# Line 119  enum ES_opgroup Line 141  enum ES_opgroup
141     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter     G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
142     G_NP1OUT,        // non-pointwise op with one output     G_NP1OUT,        // non-pointwise op with one output
143     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter     G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
144     G_TENSORPROD     // general tensor product     G_TENSORPROD,    // general tensor product
145       G_NP1OUT_2P,     // non-pointwise op with one output requiring two params
146       G_REDUCTION,     // non-pointwise unary op with a scalar output
147       G_CONDEVAL
148  };  };
149    
150    
# Line 133  string ES_opstrings[]={"UNKNOWN","IDENTI Line 158  string ES_opstrings[]={"UNKNOWN","IDENTI
158              "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",
159              "symmetric","nonsymmetric",              "symmetric","nonsymmetric",
160              "prod",              "prod",
161              "transpose", "trace"};              "transpose", "trace",
162  int ES_opcount=40;              "swapaxes",
163                "minval", "maxval",
164                "condEval"};
165    int ES_opcount=44;
166  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,
167              G_UNARY,G_UNARY,G_UNARY, //10              G_UNARY,G_UNARY,G_UNARY, //10
168              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 143  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT Line 171  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENT
171              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
172              G_NP1OUT,G_NP1OUT,              G_NP1OUT,G_NP1OUT,
173              G_TENSORPROD,              G_TENSORPROD,
174              G_NP1OUT_P, G_NP1OUT_P};              G_NP1OUT_P, G_NP1OUT_P,
175                G_NP1OUT_2P,
176                G_REDUCTION, G_REDUCTION,
177                G_CONDEVAL};
178  inline  inline
179  ES_opgroup  ES_opgroup
180  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 164  resultFS(DataAbstract_ptr left, DataAbst Line 195  resultFS(DataAbstract_ptr left, DataAbst
195    FunctionSpace r=right->getFunctionSpace();    FunctionSpace r=right->getFunctionSpace();
196    if (l!=r)    if (l!=r)
197    {    {
198      if (r.probeInterpolation(l))      signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
199        if (res==1)
200      {      {
201      return l;      return l;
202      }      }
203      if (l.probeInterpolation(r))      if (res==-1)
204      {      {
205      return r;      return r;
206      }      }
# Line 188  resultShape(DataAbstract_ptr left, DataA Line 220  resultShape(DataAbstract_ptr left, DataA
220        {        {
221          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.");
222        }        }
223    
224        if (left->getRank()==0)   // we need to allow scalar * anything        if (left->getRank()==0)   // we need to allow scalar * anything
225        {        {
226          return right->getShape();          return right->getShape();
# Line 215  resultShape(DataAbstract_ptr left, ES_op Line 248  resultShape(DataAbstract_ptr left, ES_op
248          int rank=left->getRank();          int rank=left->getRank();
249          if (axis_offset<0 || axis_offset>rank)          if (axis_offset<0 || axis_offset>rank)
250          {          {
251                 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);              stringstream e;
252              }              e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
253              for (int i=0; i<rank; i++)              throw DataException(e.str());
254            }
255            for (int i=0; i<rank; i++)
256          {          {
257             int index = (axis_offset+i)%rank;             int index = (axis_offset+i)%rank;
258                 sh.push_back(s[index]); // Append to new shape             sh.push_back(s[index]); // Append to new shape
259              }          }
260          return sh;          return sh;
261         }         }
262      break;      break;
# Line 285  resultShape(DataAbstract_ptr left, ES_op Line 320  resultShape(DataAbstract_ptr left, ES_op
320      }      }
321  }  }
322    
323    DataTypes::ShapeType
324    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
325    {
326         // This code taken from the Data.cpp swapaxes() method
327         // Some of the checks are probably redundant here
328         int axis0_tmp,axis1_tmp;
329         const DataTypes::ShapeType& s=left->getShape();
330         DataTypes::ShapeType out_shape;
331         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
332         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
333         int rank=left->getRank();
334         if (rank<2) {
335            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
336         }
337         if (axis0<0 || axis0>rank-1) {
338            stringstream e;
339            e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
340            throw DataException(e.str());
341         }
342         if (axis1<0 || axis1>rank-1) {
343            stringstream e;
344            e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
345            throw DataException(e.str());
346         }
347         if (axis0 == axis1) {
348             throw DataException("Error - Data::swapaxes: axis indices must be different.");
349         }
350         if (axis0 > axis1) {
351             axis0_tmp=axis1;
352             axis1_tmp=axis0;
353         } else {
354             axis0_tmp=axis0;
355             axis1_tmp=axis1;
356         }
357         for (int i=0; i<rank; i++) {
358           if (i == axis0_tmp) {
359              out_shape.push_back(s[axis1_tmp]);
360           } else if (i == axis1_tmp) {
361              out_shape.push_back(s[axis0_tmp]);
362           } else {
363              out_shape.push_back(s[i]);
364           }
365         }
366        return out_shape;
367    }
368    
369    
370  // determine the output shape for the general tensor product operation  // determine the output shape for the general tensor product operation
371  // the additional parameters return information required later for the product  // the additional parameters return information required later for the product
372  // 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 349  GTPShape(DataAbstract_ptr left, DataAbst Line 431  GTPShape(DataAbstract_ptr left, DataAbst
431    return shape2;    return shape2;
432  }  }
433    
 // 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  
 // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.  
 // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined  
 // multiple times  
 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 1+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)+".");  
    }  
 }  
   
   
434  }   // end anonymous namespace  }   // end anonymous namespace
435    
436    
# Line 388  opToString(ES_optype op) Line 446  opToString(ES_optype op)
446    return ES_opstrings[op];    return ES_opstrings[op];
447  }  }
448    
449    void DataLazy::LazyNodeSetup()
450    {
451    #ifdef _OPENMP
452        int numthreads=omp_get_max_threads();
453        m_samples.resize(numthreads*m_samplesize);
454        m_sampleids=new int[numthreads];
455        for (int i=0;i<numthreads;++i)
456        {
457            m_sampleids[i]=-1;  
458        }
459    #else
460        m_samples.resize(m_samplesize);
461        m_sampleids=new int[1];
462        m_sampleids[0]=-1;
463    #endif  // _OPENMP
464    }
465    
466    
467    // Creates an identity node
468  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
469      : parent(p->getFunctionSpace(),p->getShape())      : parent(p->getFunctionSpace(),p->getShape())
470        ,m_sampleids(0),
471        m_samples(1)
472  {  {
473     if (p->isLazy())     if (p->isLazy())
474     {     {
# Line 409  LAZYDEBUG(cout << "Wrapping " << dr.get( Line 487  LAZYDEBUG(cout << "Wrapping " << dr.get(
487  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
488  }  }
489    
   
   
   
490  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
491      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
492      m_op(op),      m_op(op),
493      m_axis_offset(0),      m_axis_offset(0),
494      m_transpose(0),      m_transpose(0),
495      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
496  {  {
497     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
498     {     {
499      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.");
500     }     }
# Line 435  DataLazy::DataLazy(DataAbstract_ptr left Line 510  DataLazy::DataLazy(DataAbstract_ptr left
510     }     }
511     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
512     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
513     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
514     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
515     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
516       LazyNodeSetup();
517     SIZELIMIT     SIZELIMIT
518  }  }
519    
# Line 506  LAZYDEBUG(cout << "Right " << right.get( Line 580  LAZYDEBUG(cout << "Right " << right.get(
580      m_readytype='C';      m_readytype='C';
581     }     }
582     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());    
    m_buffsRequired=calcBuffs(m_left, m_right,m_op);  
583     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
584     m_height=max(m_left->m_height,m_right->m_height)+1;     m_height=max(m_left->m_height,m_right->m_height)+1;
585       LazyNodeSetup();
586     SIZELIMIT     SIZELIMIT
587  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
588  }  }
# Line 573  DataLazy::DataLazy(DataAbstract_ptr left Line 646  DataLazy::DataLazy(DataAbstract_ptr left
646      m_readytype='C';      m_readytype='C';
647     }     }
648     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());    
    m_buffsRequired=calcBuffs(m_left, m_right,m_op);  
649     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
650     m_height=max(m_left->m_height,m_right->m_height)+1;     m_height=max(m_left->m_height,m_right->m_height)+1;
651       LazyNodeSetup();
652     SIZELIMIT     SIZELIMIT
653  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
654  }  }
# Line 604  DataLazy::DataLazy(DataAbstract_ptr left Line 676  DataLazy::DataLazy(DataAbstract_ptr left
676     }     }
677     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
678     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
679     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
680     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
681     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
682       LazyNodeSetup();
683     SIZELIMIT     SIZELIMIT
684  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
685  }  }
# Line 635  DataLazy::DataLazy(DataAbstract_ptr left Line 706  DataLazy::DataLazy(DataAbstract_ptr left
706     }     }
707     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
708     m_left=lleft;     m_left=lleft;
    m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point  
709     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());  
710     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
711     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
712       LazyNodeSetup();
713     SIZELIMIT     SIZELIMIT
714  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
715  }  }
716    
717  DataLazy::~DataLazy()  
718    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
719        : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
720        m_op(op),
721        m_axis_offset(axis0),
722        m_transpose(axis1),
723        m_tol(0)
724  {  {
725       if ((getOpgroup(op)!=G_NP1OUT_2P))
726       {
727        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
728       }
729       DataLazy_ptr lleft;
730       if (!left->isLazy())
731       {
732        lleft=DataLazy_ptr(new DataLazy(left));
733       }
734       else
735       {
736        lleft=dynamic_pointer_cast<DataLazy>(left);
737       }
738       m_readytype=lleft->m_readytype;
739       m_left=lleft;
740       m_samplesize=getNumDPPSample()*getNoValues();
741       m_children=m_left->m_children+1;
742       m_height=m_left->m_height+1;
743       LazyNodeSetup();
744       SIZELIMIT
745    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
746  }  }
747    
748    
749  int  namespace
 DataLazy::getBuffsRequired() const  
750  {  {
     return m_buffsRequired;  
 }  
751    
752        inline int max3(int a, int b, int c)
753        {
754        int t=(a>b?a:b);
755        return (t>c?t:c);
756    
757        }
758    }
759    
760  size_t  DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
761  DataLazy::getMaxSampleSize() const      : parent(left->getFunctionSpace(), left->getShape()),
762        m_op(CONDEVAL),
763        m_axis_offset(0),
764        m_transpose(0),
765        m_tol(0)
766  {  {
767      return m_maxsamplesize;  
768       DataLazy_ptr lmask;
769       DataLazy_ptr lleft;
770       DataLazy_ptr lright;
771       if (!mask->isLazy())
772       {
773        lmask=DataLazy_ptr(new DataLazy(mask));
774       }
775       else
776       {
777        lmask=dynamic_pointer_cast<DataLazy>(mask);
778       }
779       if (!left->isLazy())
780       {
781        lleft=DataLazy_ptr(new DataLazy(left));
782       }
783       else
784       {
785        lleft=dynamic_pointer_cast<DataLazy>(left);
786       }
787       if (!right->isLazy())
788       {
789        lright=DataLazy_ptr(new DataLazy(right));
790       }
791       else
792       {
793        lright=dynamic_pointer_cast<DataLazy>(right);
794       }
795       m_readytype=lmask->m_readytype;
796       if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
797       {
798        throw DataException("Programmer Error - condEval arguments must have the same readytype");
799       }
800       m_left=lleft;
801       m_right=lright;
802       m_mask=lmask;
803       m_samplesize=getNumDPPSample()*getNoValues();
804       m_children=m_left->m_children+m_right->m_children+m_mask->m_children+1;
805       m_height=max3(m_left->m_height,m_right->m_height,m_mask->m_height)+1;
806       LazyNodeSetup();
807       SIZELIMIT
808    LAZYDEBUG(cout << "(8)Lazy created with " << m_samplesize << endl;)
809  }  }
810    
811    
812    
813  size_t  DataLazy::~DataLazy()
 DataLazy::getSampleBufferSize() const  
814  {  {
815      return m_maxsamplesize*(max(1,m_buffsRequired));     delete[] m_sampleids;
816  }  }
817    
818    
819  /*  /*
820    \brief Evaluates the expression using methods on Data.    \brief Evaluates the expression using methods on Data.
821    This does the work for the collapse method.    This does the work for the collapse method.
# Line 809  DataLazy::collapseToReady() Line 955  DataLazy::collapseToReady()
955      case TRACE:      case TRACE:
956      result=left.trace(m_axis_offset);      result=left.trace(m_axis_offset);
957      break;      break;
958        case SWAP:
959        result=left.swapaxes(m_axis_offset, m_transpose);
960        break;
961        case MINVAL:
962        result=left.minval();
963        break;
964        case MAXVAL:
965        result=left.minval();
966        break;
967      default:      default:
968      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)+".");
969    }    }
# Line 836  DataLazy::collapse() Line 991  DataLazy::collapse()
991    m_op=IDENTITY;    m_op=IDENTITY;
992  }  }
993    
994  /*  
995    \brief Compute the value of the expression (unary operation) for the given sample.  
996    \return Vector which stores the value of the subexpression for the given sample.  
997    \param v A vector to store intermediate results.  
998    \param offset Index in v to begin storing results.  
999    \param sampleNo Sample number to evaluate.  #define PROC_OP(TYPE,X)                               \
1000    \param roffset (output parameter) the offset in the return vector where the result begins.      for (int j=0;j<onumsteps;++j)\
1001        {\
1002    The return value will be an existing vector so do not deallocate it.        for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1003    If the result is stored in v it should be stored at the offset given.        { \
1004    Everything from offset to the end of v should be considered available for this method to use.  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1005  */  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1006  DataTypes::ValueType*           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1007  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1008             lroffset+=leftstep; \
1009             rroffset+=rightstep; \
1010          }\
1011          lroffset+=oleftstep;\
1012          rroffset+=orightstep;\
1013        }
1014    
1015    
1016    // The result will be stored in m_samples
1017    // The return value is a pointer to the DataVector, offset is the offset within the return value
1018    const DataTypes::ValueType*
1019    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
1020    {
1021    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1022        // collapse so we have a 'E' node or an IDENTITY for some other type
1023      if (m_readytype!='E' && m_op!=IDENTITY)
1024      {
1025        collapse();
1026      }
1027      if (m_op==IDENTITY)  
1028      {
1029        const ValueType& vec=m_id->getVectorRO();
1030        roffset=m_id->getPointOffset(sampleNo, 0);
1031    #ifdef LAZY_STACK_PROF
1032    int x;
1033    if (&x<stackend[omp_get_thread_num()])
1034    {
1035           stackend[omp_get_thread_num()]=&x;
1036    }
1037    #endif
1038        return &(vec);
1039      }
1040      if (m_readytype!='E')
1041      {
1042        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1043      }
1044      if (m_sampleids[tid]==sampleNo)
1045      {
1046        roffset=tid*m_samplesize;
1047        return &(m_samples);        // sample is already resolved
1048      }
1049      m_sampleids[tid]=sampleNo;
1050    
1051      switch (getOpgroup(m_op))
1052      {
1053      case G_UNARY:
1054      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1055      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1056      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1057      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1058      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1059      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1060      case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
1061      case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
1062      default:
1063        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1064      }
1065    }
1066    
1067    const DataTypes::ValueType*
1068    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
1069  {  {
1070      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1071      // 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
1072      // processing single points.      // processing single points.
1073        // we will also know we won't get identity nodes
1074    if (m_readytype!='E')    if (m_readytype!='E')
1075    {    {
1076      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1077    }    }
1078    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);    if (m_op==IDENTITY)
1079    const double* left=&((*vleft)[roffset]);    {
1080    double* result=&(v[offset]);      throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1081    roffset=offset;    }
1082      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1083      const double* left=&((*leftres)[roffset]);
1084      roffset=m_samplesize*tid;
1085      double* result=&(m_samples[roffset]);
1086    switch (m_op)    switch (m_op)
1087    {    {
1088      case SIN:        case SIN:  
# Line 971  DataLazy::resolveUnary(ValueType& v, siz Line 1192  DataLazy::resolveUnary(ValueType& v, siz
1192      default:      default:
1193      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1194    }    }
1195    return &v;    return &(m_samples);
1196  }  }
1197    
1198    
1199    const DataTypes::ValueType*
1200    DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)
1201    {
1202        // we assume that any collapsing has been done before we get here
1203        // since we only have one argument we don't need to think about only
1204        // processing single points.
1205        // we will also know we won't get identity nodes
1206      if (m_readytype!='E')
1207      {
1208        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1209      }
1210      if (m_op==IDENTITY)
1211      {
1212        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1213      }
1214      size_t loffset=0;
1215      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1216    
1217      roffset=m_samplesize*tid;
1218      unsigned int ndpps=getNumDPPSample();
1219      unsigned int psize=DataTypes::noValues(m_left->getShape());
1220      double* result=&(m_samples[roffset]);
1221      switch (m_op)
1222      {
1223        case MINVAL:
1224        {
1225          for (unsigned int z=0;z<ndpps;++z)
1226          {
1227            FMin op;
1228            *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1229            loffset+=psize;
1230            result++;
1231          }
1232        }
1233        break;
1234        case MAXVAL:
1235        {
1236          for (unsigned int z=0;z<ndpps;++z)
1237          {
1238          FMax op;
1239          *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1240          loffset+=psize;
1241          result++;
1242          }
1243        }
1244        break;
1245        default:
1246        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1247      }
1248      return &(m_samples);
1249    }
1250    
1251    const DataTypes::ValueType*
1252    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  
1253  {  {
1254      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1255      // 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
1256      // processing single points.      // processing single points.
1257    if (m_readytype!='E')    if (m_readytype!='E')
1258    {    {
1259      throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");      throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1260    }    }
1261      // since we can't write the result over the input, we need a result offset further along    if (m_op==IDENTITY)
1262    size_t subroffset=roffset+m_samplesize;    {
1263  LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1264    const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);    }
1265    roffset=offset;    size_t subroffset;
1266      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1267      roffset=m_samplesize*tid;
1268    size_t loop=0;    size_t loop=0;
1269    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1270    size_t step=getNoValues();    size_t step=getNoValues();
1271      size_t offset=roffset;
1272    switch (m_op)    switch (m_op)
1273    {    {
1274      case SYM:      case SYM:
1275      for (loop=0;loop<numsteps;++loop)      for (loop=0;loop<numsteps;++loop)
1276      {      {
1277          DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);          DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1278          subroffset+=step;          subroffset+=step;
1279          offset+=step;          offset+=step;
1280      }      }
# Line 1022  LAZYDEBUG(cerr << "subroffset=" << subro Line 1282  LAZYDEBUG(cerr << "subroffset=" << subro
1282      case NSYM:      case NSYM:
1283      for (loop=0;loop<numsteps;++loop)      for (loop=0;loop<numsteps;++loop)
1284      {      {
1285          DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);          DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1286          subroffset+=step;          subroffset+=step;
1287          offset+=step;          offset+=step;
1288      }      }
# Line 1030  LAZYDEBUG(cerr << "subroffset=" << subro Line 1290  LAZYDEBUG(cerr << "subroffset=" << subro
1290      default:      default:
1291      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1292    }    }
1293    return &v;    return &m_samples;
1294  }  }
1295    
1296  /*  const DataTypes::ValueType*
1297    \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  
1298  {  {
1299      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
1300      // 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
1301      // processing single points.      // processing single points.
1302    if (m_readytype!='E')    if (m_readytype!='E')
1303    {    {
1304      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.");
1305      }
1306      if (m_op==IDENTITY)
1307      {
1308        throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1309    }    }
     // since we can't write the result over the input, we need a result offset further along  
1310    size_t subroffset;    size_t subroffset;
1311    const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);    size_t offset;
1312  LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1313  LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)    roffset=m_samplesize*tid;
1314    roffset=offset;    offset=roffset;
1315    size_t loop=0;    size_t loop=0;
1316    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1317    size_t outstep=getNoValues();    size_t outstep=getNoValues();
1318    size_t instep=m_left->getNoValues();    size_t instep=m_left->getNoValues();
 LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)  
1319    switch (m_op)    switch (m_op)
1320    {    {
1321      case TRACE:      case TRACE:
1322      for (loop=0;loop<numsteps;++loop)      for (loop=0;loop<numsteps;++loop)
1323      {      {
1324  size_t zz=sampleNo*getNumDPPSample()+loop;              DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
 if (zz==5800)  
 {  
 LAZYDEBUG(cerr << "point=" <<  zz<< endl;)  
 LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)  
 LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)  
 LAZYDEBUG(cerr << subroffset << endl;)  
 LAZYDEBUG(cerr << "output=" << offset << endl;)  
 }  
             DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);  
 if (zz==5800)  
 {  
 LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)  
 }  
1325          subroffset+=instep;          subroffset+=instep;
1326          offset+=outstep;          offset+=outstep;
1327      }      }
# Line 1092  LAZYDEBUG(cerr << "Result of trace=" << Line 1329  LAZYDEBUG(cerr << "Result of trace=" <<
1329      case TRANS:      case TRANS:
1330      for (loop=0;loop<numsteps;++loop)      for (loop=0;loop<numsteps;++loop)
1331      {      {
1332              DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);              DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1333          subroffset+=instep;          subroffset+=instep;
1334          offset+=outstep;          offset+=outstep;
1335      }      }
# Line 1100  LAZYDEBUG(cerr << "Result of trace=" << Line 1337  LAZYDEBUG(cerr << "Result of trace=" <<
1337      default:      default:
1338      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1339    }    }
1340    return &v;    return &m_samples;
1341  }  }
1342    
1343    
1344  #define PROC_OP(TYPE,X)                               \  const DataTypes::ValueType*
1345      for (int j=0;j<onumsteps;++j)\  DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1346      {\  {
1347        for (int i=0;i<numsteps;++i,resultp+=resultStep) \    if (m_readytype!='E')
1348        { \    {
1349  LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\      throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1350  LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\    }
1351           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \    if (m_op==IDENTITY)
1352  LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \    {
1353           lroffset+=leftstep; \      throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1354           rroffset+=rightstep; \    }
1355        }\    size_t subroffset;
1356        lroffset+=oleftstep;\    size_t offset;
1357        rroffset+=orightstep;\    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1358      roffset=m_samplesize*tid;
1359      offset=roffset;
1360      size_t loop=0;
1361      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1362      size_t outstep=getNoValues();
1363      size_t instep=m_left->getNoValues();
1364      switch (m_op)
1365      {
1366        case SWAP:
1367        for (loop=0;loop<numsteps;++loop)
1368        {
1369                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1370            subroffset+=instep;
1371            offset+=outstep;
1372      }      }
1373        break;
1374        default:
1375        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1376      }
1377      return &m_samples;
1378    }
1379    
1380    const DataTypes::ValueType*
1381    DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset)
1382    {
1383      if (m_readytype!='E')
1384      {
1385        throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1386      }
1387      if (m_op!=CONDEVAL)
1388      {
1389        throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1390      }
1391      size_t subroffset;
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      for (int i=0;i<m_samplesize;++i)
1408      {
1409        m_samples[roffset+i]=(*srcres)[subroffset+i];  
1410      }
1411    
1412      return &m_samples;
1413    }
1414    
 /*  
   \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.  
 */  
1415  // 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
1416  // 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.
1417  // 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 1141  LAZYDEBUG(cout << " result=      " << re Line 1421  LAZYDEBUG(cout << " result=      " << re
1421  // There is an additional complication when scalar operations are considered.  // There is an additional complication when scalar operations are considered.
1422  // For example, 2+Vector.  // For example, 2+Vector.
1423  // In this case each double within the point is treated individually  // In this case each double within the point is treated individually
1424  DataTypes::ValueType*  const DataTypes::ValueType*
1425  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1426  {  {
1427  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)  LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1428    
# Line 1254  LAZYDEBUG(cout << "Resolve binary: " << Line 1534  LAZYDEBUG(cout << "Resolve binary: " <<
1534    
1535    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1536      // Get the values of sub-expressions      // Get the values of sub-expressions
1537    const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
1538      // calcBufss for why we can't put offset as the 2nd param above    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
   const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note  
     // the right child starts further along.  
1539  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1540  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1541  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
# Line 1266  LAZYDEBUG(cout << "onumsteps=" << onumst Line 1544  LAZYDEBUG(cout << "onumsteps=" << onumst
1544  LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)  LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1545  LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)  LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1546    
1547    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1548    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1549    
1550    
1551    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    roffset=m_samplesize*tid;
1552      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
1553    switch(m_op)    switch(m_op)
1554    {    {
1555      case ADD:      case ADD:
# Line 1288  LAZYDEBUG(cout << "" << LS << RS << LN < Line 1570  LAZYDEBUG(cout << "" << LS << RS << LN <
1570      default:      default:
1571      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1572    }    }
1573    roffset=offset;  LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1574    return &v;    return &m_samples;
1575  }  }
1576    
1577    
   
 /*  
   \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.  
 */  
1578  // 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
1579  // 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.
1580  // unlike the other resolve helpers, we must treat these datapoints separately.  // unlike the other resolve helpers, we must treat these datapoints separately.
1581  DataTypes::ValueType*  const DataTypes::ValueType*
1582  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
1583  {  {
1584  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1585    
1586    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
1587      // first work out which of the children are expanded      // first work out which of the children are expanded
1588    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1589    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1590    int steps=getNumDPPSample();    int steps=getNumDPPSample();
 /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);  
   int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/  
1591    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method    int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1592    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1593    
1594    int resultStep=getNoValues();    int resultStep=getNoValues();
1595      // Get the values of sub-expressions (leave a gap of one sample for the result).    roffset=m_samplesize*tid;
1596    int gap=offset+m_samplesize;      size_t offset=roffset;
   
 LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)  
   
   const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);  
   gap+=m_left->getMaxSampleSize();  
   
1597    
1598  LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)    const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1599    
1600      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
   const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);  
1601    
1602  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1603  cout << getNoValues() << endl;)  cout << getNoValues() << endl;)
 LAZYDEBUG(cerr << "Result of left=";)  
 LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;  
1604    
1605  for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)  
 {  
 cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";  
 if (i%4==0) cout << endl;  
 })  
 LAZYDEBUG(cerr << "\nResult of right=" << endl;)  
 LAZYDEBUG(  
 for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)  
 {  
 cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";  
 if (i%4==0) cout << endl;  
 }  
 cerr << endl;  
 )  
1606  LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)  LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1607  LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)  LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1608  LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)  LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
# Line 1366  LAZYDEBUG(cout << "m_samplesize=" << m_s Line 1611  LAZYDEBUG(cout << "m_samplesize=" << m_s
1611  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1612  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)  LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1613    
1614    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
1615    switch(m_op)    switch(m_op)
1616    {    {
1617      case PROD:      case PROD:
1618      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1619      {      {
   
 LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)  
 LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)  
 LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)  
   
1620            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1621            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1622    
# Line 1385  LAZYDEBUG(cout << DataTypes::pointToStri Line 1625  LAZYDEBUG(cout << DataTypes::pointToStri
1625    
1626            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);
1627    
 LAZYDEBUG(cout << "Results--\n";  
 {  
   DataVector dv(getNoValues());  
 for (int z=0;z<getNoValues();++z)  
 {  
   cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";  
   if (z%4==0) cout << endl;  
   dv[z]=resultp[z];  
 }  
 cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");  
 cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;  
 }  
 )  
1628        lroffset+=leftStep;        lroffset+=leftStep;
1629        rroffset+=rightStep;        rroffset+=rightStep;
1630      }      }
# Line 1406  cout << "\nWritten to: " << resultp << " Line 1633  cout << "\nWritten to: " << resultp << "
1633      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1634    }    }
1635    roffset=offset;    roffset=offset;
1636    return &v;    return &m_samples;
1637  }  }
1638    
1639    
   
 /*  
   \brief Compute the value of the expression for the given sample.  
   \return Vector which stores the value of the subexpression for the given sample.  
   \param v A vector to store intermediate results.  
   \param offset Index in v to begin storing results.  
   \param sampleNo Sample number to evaluate.  
   \param roffset (output parameter) the offset in the return vector where the result begins.  
   
   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.  
   
 // the roffset is the offset within the returned vector where the data begins  
 const DataTypes::ValueType*  
 DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  
 {  
 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)  
     // collapse so we have a 'E' node or an IDENTITY for some other type  
   if (m_readytype!='E' && m_op!=IDENTITY)  
   {  
     collapse();  
   }  
   if (m_op==IDENTITY)    
   {  
     const ValueType& vec=m_id->getVectorRO();  
     if (m_readytype=='C')  
     {  
     roffset=0;  
 LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)  
     return &(vec);  
     }  
     roffset=m_id->getPointOffset(sampleNo, 0);  
 LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)  
     return &(vec);  
   }  
   if (m_readytype!='E')  
   {  
     throw DataException("Programmer Error - Collapse did not produce an expanded node.");  
   }  
   switch (getOpgroup(m_op))  
   {  
   case G_UNARY:  
   case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);  
   case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);  
   case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);  
   case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);  
   case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);  
   default:  
     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
   }  
   
 }  
   
1640  const DataTypes::ValueType*  const DataTypes::ValueType*
1641  DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)  DataLazy::resolveSample(int sampleNo, size_t& roffset)
1642  {  {
1643  #ifdef _OPENMP  #ifdef _OPENMP
1644      int tid=omp_get_thread_num();      int tid=omp_get_thread_num();
1645  #else  #else
1646      int tid=0;      int tid=0;
1647  #endif  #endif
1648      return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);  
1649    #ifdef LAZY_STACK_PROF
1650        stackstart[tid]=&tid;
1651        stackend[tid]=&tid;
1652        const DataTypes::ValueType* r=resolveNodeSample(tid, sampleNo, roffset);
1653        size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1654        #pragma omp critical
1655        if (d>maxstackuse)
1656        {
1657    cout << "Max resolve Stack use " << d << endl;
1658            maxstackuse=d;
1659        }
1660        return r;
1661    #else
1662        return resolveNodeSample(tid, sampleNo, roffset);
1663    #endif
1664  }  }
1665    
1666    
1667  // This needs to do the work of the idenity constructor  // This needs to do the work of the identity constructor
1668  void  void
1669  DataLazy::resolveToIdentity()  DataLazy::resolveToIdentity()
1670  {  {
1671     if (m_op==IDENTITY)     if (m_op==IDENTITY)
1672      return;      return;
1673     DataReady_ptr p=resolve();     DataReady_ptr p=resolveNodeWorker();
1674     makeIdentity(p);     makeIdentity(p);
1675  }  }
1676    
# Line 1501  void DataLazy::makeIdentity(const DataRe Line 1686  void DataLazy::makeIdentity(const DataRe
1686     else if(p->isExpanded()) {m_readytype='E';}     else if(p->isExpanded()) {m_readytype='E';}
1687     else if (p->isTagged()) {m_readytype='T';}     else if (p->isTagged()) {m_readytype='T';}
1688     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
    m_buffsRequired=1;  
1689     m_samplesize=p->getNumDPPSample()*p->getNoValues();     m_samplesize=p->getNumDPPSample()*p->getNoValues();
    m_maxsamplesize=m_samplesize;  
1690     m_left.reset();     m_left.reset();
1691     m_right.reset();     m_right.reset();
1692  }  }
1693    
1694  // To simplify the memory management, all threads operate on one large vector, rather than one each.  
 // Each sample is evaluated independently and copied into the result DataExpanded.  
1695  DataReady_ptr  DataReady_ptr
1696  DataLazy::resolve()  DataLazy::resolve()
1697  {  {
1698        resolveToIdentity();
1699        return m_id;
1700    }
1701    
1702    
1703    /* This is really a static method but I think that caused problems in windows */
1704    void
1705    DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
1706    {
1707      if (dats.empty())
1708      {
1709        return;
1710      }
1711      vector<DataLazy*> work;
1712      FunctionSpace fs=dats[0]->getFunctionSpace();
1713      bool match=true;
1714      for (int i=dats.size()-1;i>=0;--i)
1715      {
1716        if (dats[i]->m_readytype!='E')
1717        {
1718            dats[i]->collapse();
1719        }
1720        if (dats[i]->m_op!=IDENTITY)
1721        {
1722            work.push_back(dats[i]);
1723            if (fs!=dats[i]->getFunctionSpace())
1724            {
1725                match=false;
1726            }
1727        }
1728      }
1729      if (work.empty())
1730      {
1731        return;     // no work to do
1732      }
1733      if (match)    // all functionspaces match.  Yes I realise this is overly strict
1734      {     // it is possible that dats[0] is one of the objects which we discarded and
1735            // all the other functionspaces match.
1736        vector<DataExpanded*> dep;
1737        vector<ValueType*> vecs;
1738        for (int i=0;i<work.size();++i)
1739        {
1740            dep.push_back(new DataExpanded(fs,work[i]->getShape(), ValueType(work[i]->getNoValues())));
1741            vecs.push_back(&(dep[i]->getVectorRW()));
1742        }
1743        int totalsamples=work[0]->getNumSamples();
1744        const ValueType* res=0; // Storage for answer
1745        int sample;
1746        #pragma omp parallel private(sample, res)
1747        {
1748            size_t roffset=0;
1749            #pragma omp for schedule(static)
1750            for (sample=0;sample<totalsamples;++sample)
1751            {
1752            roffset=0;
1753            int j;
1754            for (j=work.size()-1;j>=0;--j)
1755            {
1756    #ifdef _OPENMP
1757                    res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1758    #else
1759                    res=work[j]->resolveNodeSample(0,sample,roffset);
1760    #endif
1761                    DataVector::size_type outoffset=dep[j]->getPointOffset(sample,0);
1762                    memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(DataVector::ElementType));
1763            }
1764            }
1765        }
1766        // Now we need to load the new results as identity ops into the lazy nodes
1767        for (int i=work.size()-1;i>=0;--i)
1768        {
1769            work[i]->makeIdentity(boost::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1770        }
1771      }
1772      else  // functionspaces do not match
1773      {
1774        for (int i=0;i<work.size();++i)
1775        {
1776            work[i]->resolveToIdentity();
1777        }
1778      }
1779    }
1780    
1781  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)  
1782  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)  
1783    // This version of resolve uses storage in each node to hold results
1784    DataReady_ptr
1785    DataLazy::resolveNodeWorker()
1786    {
1787    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
1788    {    {
1789      collapse();      collapse();
# Line 1525  LAZYDEBUG(cout << "Buffers=" << m_buffsR Line 1793  LAZYDEBUG(cout << "Buffers=" << m_buffsR
1793      return m_id;      return m_id;
1794    }    }
1795      // 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=omp_get_max_threads();  
 #endif  
   ValueType v(numthreads*threadbuffersize);  
 LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)  
1796    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1797    ValueType& resvec=result->getVectorRW();    ValueType& resvec=result->getVectorRW();
1798    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1799    
1800    int sample;    int sample;
   size_t outoffset;     // offset in the output data  
1801    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1802    const ValueType* res=0;   // Vector storing the answer    const ValueType* res=0;   // Storage for answer
   size_t resoffset=0;       // where in the vector to find the answer  
1803  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1804    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel private(sample,res)
   for (sample=0;sample<totalsamples;++sample)  
1805    {    {
1806        if (sample==0)  {ENABLEDEBUG}      size_t roffset=0;
1807    #ifdef LAZY_STACK_PROF
1808  //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}      stackstart[omp_get_thread_num()]=&roffset;
1809  LAZYDEBUG(cout << "################################# " << sample << endl;)      stackend[omp_get_thread_num()]=&roffset;
1810    #endif
1811        #pragma omp for schedule(static)
1812        for (sample=0;sample<totalsamples;++sample)
1813        {
1814            roffset=0;
1815  #ifdef _OPENMP  #ifdef _OPENMP
1816      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);              res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1817  #else  #else
1818      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.              res=resolveNodeSample(0,sample,roffset);
1819  #endif  #endif
1820  LAZYDEBUG(cerr << "-------------------------------- " << endl;)  LAZYDEBUG(cout << "Sample #" << sample << endl;)
1821  LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)  LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1822      outoffset=result->getPointOffset(sample,0);              DataVector::size_type outoffset=result->getPointOffset(sample,0);
1823  LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)              memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
1824      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector      }
1825      {    }
1826  LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)  #ifdef LAZY_STACK_PROF
1827      resvec[outoffset]=(*res)[resoffset];    for (int i=0;i<getNumberOfThreads();++i)
1828      }    {
1829  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)      size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1830  LAZYDEBUG(cerr << "*********************************" << endl;)  //  cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " <<  r << endl;
1831      DISABLEDEBUG      if (r>maxstackuse)
1832        {
1833            maxstackuse=r;
1834        }
1835    }    }
1836      cout << "Max resolve Stack use=" << maxstackuse << endl;
1837    #endif
1838    return resptr;    return resptr;
1839  }  }
1840    
# Line 1574  std::string Line 1842  std::string
1842  DataLazy::toString() const  DataLazy::toString() const
1843  {  {
1844    ostringstream oss;    ostringstream oss;
1845    oss << "Lazy Data:";    oss << "Lazy Data: [depth=" << m_height<< "] ";
1846    intoString(oss);    switch (escriptParams.getLAZY_STR_FMT())
1847      {
1848      case 1:   // tree format
1849        oss << endl;
1850        intoTreeString(oss,"");
1851        break;
1852      case 2:   // just the depth
1853        break;
1854      default:
1855        intoString(oss);
1856        break;
1857      }
1858    return oss.str();    return oss.str();
1859  }  }
1860    
# Line 1616  DataLazy::intoString(ostringstream& oss) Line 1895  DataLazy::intoString(ostringstream& oss)
1895    case G_UNARY_P:    case G_UNARY_P:
1896    case G_NP1OUT:    case G_NP1OUT:
1897    case G_NP1OUT_P:    case G_NP1OUT_P:
1898      case G_REDUCTION:
1899      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1900      m_left->intoString(oss);      m_left->intoString(oss);
1901      oss << ')';      oss << ')';
# Line 1627  DataLazy::intoString(ostringstream& oss) Line 1907  DataLazy::intoString(ostringstream& oss)
1907      m_right->intoString(oss);      m_right->intoString(oss);
1908      oss << ')';      oss << ')';
1909      break;      break;
1910      case G_NP1OUT_2P:
1911        oss << opToString(m_op) << '(';
1912        m_left->intoString(oss);
1913        oss << ", " << m_axis_offset << ", " << m_transpose;
1914        oss << ')';
1915        break;
1916      case G_CONDEVAL:
1917        oss << opToString(m_op)<< '(' ;
1918        m_mask->intoString(oss);
1919        oss << " ? ";
1920        m_left->intoString(oss);
1921        oss << " : ";
1922        m_right->intoString(oss);
1923        oss << ')';
1924        break;
1925      default:
1926        oss << "UNKNOWN";
1927      }
1928    }
1929    
1930    
1931    void
1932    DataLazy::intoTreeString(ostringstream& oss, string indent) const
1933    {
1934      oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1935      switch (getOpgroup(m_op))
1936      {
1937      case G_IDENTITY:
1938        if (m_id->isExpanded())
1939        {
1940           oss << "E";
1941        }
1942        else if (m_id->isTagged())
1943        {
1944          oss << "T";
1945        }
1946        else if (m_id->isConstant())
1947        {
1948          oss << "C";
1949        }
1950        else
1951        {
1952          oss << "?";
1953        }
1954        oss << '@' << m_id.get() << endl;
1955        break;
1956      case G_BINARY:
1957        oss << opToString(m_op) << endl;
1958        indent+='.';
1959        m_left->intoTreeString(oss, indent);
1960        m_right->intoTreeString(oss, indent);
1961        break;
1962      case G_UNARY:
1963      case G_UNARY_P:
1964      case G_NP1OUT:
1965      case G_NP1OUT_P:
1966      case G_REDUCTION:
1967        oss << opToString(m_op) << endl;
1968        indent+='.';
1969        m_left->intoTreeString(oss, indent);
1970        break;
1971      case G_TENSORPROD:
1972        oss << opToString(m_op) << endl;
1973        indent+='.';
1974        m_left->intoTreeString(oss, indent);
1975        m_right->intoTreeString(oss, indent);
1976        break;
1977      case G_NP1OUT_2P:
1978        oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1979        indent+='.';
1980        m_left->intoTreeString(oss, indent);
1981        break;
1982    default:    default:
1983      oss << "UNKNOWN";      oss << "UNKNOWN";
1984    }    }
1985  }  }
1986    
1987    
1988  DataAbstract*  DataAbstract*
1989  DataLazy::deepCopy()  DataLazy::deepCopy()
1990  {  {
1991    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1992    {    {
1993    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1994    case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);    case G_UNARY:
1995      case G_REDUCTION:      return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1996      case G_UNARY_P:   return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1997    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);
1998    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);
1999    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);
2000      case G_NP1OUT_P:   return new DataLazy(m_left->deepCopy()->getPtr(),m_op,  m_axis_offset);
2001      case G_NP1OUT_2P:  return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2002    default:    default:
2003      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)+".");
2004    }    }
2005  }  }
2006    
2007    
2008    
2009  // There is no single, natural interpretation of getLength on DataLazy.  // There is no single, natural interpretation of getLength on DataLazy.
2010  // Instances of DataReady can look at the size of their vectors.  // Instances of DataReady can look at the size of their vectors.
2011  // 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 1736  DataLazy::setToZero() Line 2094  DataLazy::setToZero()
2094  //   m_readytype='C';  //   m_readytype='C';
2095  //   m_buffsRequired=1;  //   m_buffsRequired=1;
2096    
2097      privdebug=privdebug;  // to stop the compiler complaining about unused privdebug
2098    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
   
2099  }  }
2100    
2101  bool  bool

Legend:
Removed from v.2472  
changed lines
  Added in v.4264

  ViewVC Help
Powered by ViewVC 1.1.26