/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataLazy.cpp

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

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

Legend:
Removed from v.2157  
changed lines
  Added in v.2500

  ViewVC Help
Powered by ViewVC 1.1.26