/[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 1868 by jfenwick, Thu Oct 9 06:30:49 2008 UTC revision 1888 by jfenwick, Wed Oct 15 04:00:21 2008 UTC
# Line 19  Line 19 
19  #ifdef PASO_MPI  #ifdef PASO_MPI
20  #include <mpi.h>  #include <mpi.h>
21  #endif  #endif
22    #ifdef _OPENMP
23    #include <omp.h>
24    #endif
25  #include "FunctionSpace.h"  #include "FunctionSpace.h"
26  #include "DataTypes.h"  #include "DataTypes.h"
27  #include "Data.h"  #include "Data.h"
28    #include "UnaryFuncs.h"     // for escript::fsign
29    
30  using namespace std;  using namespace std;
31  using namespace boost;  using namespace boost;
# Line 35  opToString(ES_optype op); Line 39  opToString(ES_optype op);
39  namespace  namespace
40  {  {
41    
42    
43    
44    enum ES_opgroup
45    {
46       G_UNKNOWN,
47       G_IDENTITY,
48       G_BINARY,
49       G_UNARY
50    };
51    
52    
53    
54    
55    string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",
56                "asin","acos","atan","sinh","cosh","tanh","erf",
57                "asinh","acosh","atanh",
58                "log10","log","sign","abs","neg","pos","exp","sqrt",
59                "1/","where>0","where<0","where>=0","where<=0"};
60    int ES_opcount=32;
61    ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9
62                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16
63                G_UNARY,G_UNARY,G_UNARY,                    // 19
64                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27
65                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
66    inline
67    ES_opgroup
68    getOpgroup(ES_optype op)
69    {
70      return opgroups[op];
71    }
72    
73  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
74  FunctionSpace  FunctionSpace
75  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
# Line 43  resultFS(DataAbstract_ptr left, DataAbst Line 78  resultFS(DataAbstract_ptr left, DataAbst
78      // maybe we need an interpolate node -      // maybe we need an interpolate node -
79      // that way, if interpolate is required in any other op we can just throw a      // that way, if interpolate is required in any other op we can just throw a
80      // programming error exception.      // programming error exception.
81    
82    
83        if (left->getFunctionSpace()!=right->getFunctionSpace())
84        {
85            throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
86        }
87      return left->getFunctionSpace();      return left->getFunctionSpace();
88  }  }
89    
# Line 50  resultFS(DataAbstract_ptr left, DataAbst Line 91  resultFS(DataAbstract_ptr left, DataAbst
91  DataTypes::ShapeType  DataTypes::ShapeType
92  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
93  {  {
94      return DataTypes::scalarShape;      if (left->getShape()!=right->getShape())
95        {
96            throw DataException("Shapes not the same - shapes must match for lazy data.");
97        }
98        return left->getShape();
99  }  }
100    
101  size_t  size_t
102  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
103  {  {
104     switch(op)     switch (getOpgroup(op))
105     {     {
106     case IDENTITY: return left->getLength();     case G_BINARY: return left->getLength();
107     case ADD:    // the length is preserved in these ops     case G_UNARY: return left->getLength();
    case SUB:  
    case MUL:  
    case DIV: return left->getLength();  
108     default:     default:
109      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
110       }
111    }
112    
113    int
114    calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
115    {
116       switch(getOpgroup(op))
117       {
118       case G_IDENTITY: return 1;
119       case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
120       case G_UNARY: return max(left->getBuffsRequired(),1);
121       default:
122        throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
123     }     }
124  }  }
125    
126  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};  
 int ES_opcount=5;  
127    
128  }   // end anonymous namespace  }   // end anonymous namespace
129    
# Line 88  opToString(ES_optype op) Line 141  opToString(ES_optype op)
141    
142  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
143      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape()),
     m_left(p),  
144      m_op(IDENTITY)      m_op(IDENTITY)
145  {  {
146     length=resultLength(m_left,m_right,m_op);     if (p->isLazy())
147       {
148        // TODO: fix this.   We could make the new node a copy of p?
149        // I don't want identity of Lazy.
150        // Question: Why would that be so bad?
151        // Answer: We assume that the child of ID is something we can call getVector on
152        throw DataException("Programmer error - attempt to create identity from a DataLazy.");
153       }
154       else
155       {
156        m_id=dynamic_pointer_cast<DataReady>(p);
157       }
158       m_length=p->getLength();
159       m_buffsRequired=1;
160       m_samplesize=getNumDPPSample()*getNoValues();
161    cout << "(1)Lazy created with " << m_samplesize << endl;
162    }
163    
164    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
165        : parent(left->getFunctionSpace(),left->getShape()),
166        m_op(op)
167    {
168       if (getOpgroup(op)!=G_UNARY)
169       {
170        throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
171       }
172       DataLazy_ptr lleft;
173       if (!left->isLazy())
174       {
175        lleft=DataLazy_ptr(new DataLazy(left));
176       }
177       else
178       {
179        lleft=dynamic_pointer_cast<DataLazy>(left);
180       }
181       m_length=left->getLength();
182       m_left=lleft;
183       m_buffsRequired=1;
184       m_samplesize=getNumDPPSample()*getNoValues();
185  }  }
186    
187  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
188    DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
189      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
190      m_left(left),      m_left(left),
191      m_right(right),      m_right(right),
192      m_op(op)      m_op(op)
193  {  {
194     length=resultLength(m_left,m_right,m_op);     if (getOpgroup(op)!=G_BINARY)
195       {
196        throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
197       }
198       m_length=resultLength(m_left,m_right,m_op);
199       m_samplesize=getNumDPPSample()*getNoValues();
200       m_buffsRequired=calcBuffs(m_left, m_right, m_op);
201    cout << "(2)Lazy created with " << m_samplesize << endl;
202    }
203    
204    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
205        : parent(resultFS(left,right,op), resultShape(left,right,op)),
206        m_op(op)
207    {
208       if (getOpgroup(op)!=G_BINARY)
209       {
210        throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");
211       }
212       if (left->isLazy())
213       {
214        m_left=dynamic_pointer_cast<DataLazy>(left);
215       }
216       else
217       {
218        m_left=DataLazy_ptr(new DataLazy(left));
219       }
220       if (right->isLazy())
221       {
222        m_right=dynamic_pointer_cast<DataLazy>(right);
223       }
224       else
225       {
226        m_right=DataLazy_ptr(new DataLazy(right));
227       }
228    
229       m_length=resultLength(m_left,m_right,m_op);
230       m_samplesize=getNumDPPSample()*getNoValues();
231       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
232    cout << "(3)Lazy created with " << m_samplesize << endl;
233  }  }
234    
235    
236  DataLazy::~DataLazy()  DataLazy::~DataLazy()
237  {  {
238  }  }
239    
240  // If resolving records a pointer to the resolved Data we may need to rethink the const on this method  
241  DataReady_ptr  int
242  DataLazy::resolve()  DataLazy::getBuffsRequired() const
243  {  {
244    DataReady_ptr left;      return m_buffsRequired;
245    DataReady_ptr right;  }
246    if (m_left.get()!=0)  
247    {  
248      left=m_left->resolve();  // the vector and the offset are a place where the method could write its data if it wishes
249    }  // it is not obligated to do so. For example, if it has its own storage already, it can use that.
250    if (m_right.get()!=0)  // Hence the return value to indicate where the data is actually stored.
251    {  // Regardless, the storage should be assumed to be used, even if it isn't.
252      right=m_right->resolve();  const double*
253    }  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const
254    switch (m_op)  {
255    {    if (m_op==IDENTITY)  
256      case IDENTITY: return left;    {
257      case ADD:      const ValueType& vec=m_id->getVector();
258      // Hmm we could get interpolation here, better be careful      return &(vec[m_id->getPointOffset(sampleNo, 0)]);
259        return C_TensorBinaryOperation(Data(left),Data(right),plus<double>()).borrowReadyPtr();    }
260      case SUB:    size_t rightoffset=offset+m_samplesize;
261        return C_TensorBinaryOperation(Data(left),Data(right),minus<double>()).borrowReadyPtr();    const double* left=m_left->resolveSample(v,sampleNo,offset);
262      case MUL:    const double* right=0;
263        return C_TensorBinaryOperation(Data(left),Data(right),multiplies<double>()).borrowReadyPtr();    if (getOpgroup(m_op)==G_BINARY)
264      case DIV:    {
265        return C_TensorBinaryOperation(Data(left),Data(right),divides<double>()).borrowReadyPtr();      right=m_right->resolveSample(v,sampleNo,rightoffset);
266      }
267      double* result=&(v[offset]);
268      {
269        switch(m_op)
270        {
271        case ADD:       // since these are pointwise ops, pretend each sample is one point
272        tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
273        break;
274        case SUB:      
275        tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
276        break;
277        case MUL:      
278        tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
279        break;
280        case DIV:      
281        tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
282        break;
283    // unary ops
284        case SIN:
285        tensor_unary_operation(m_samplesize, left, result, ::sin);
286        break;
287        case COS:
288        tensor_unary_operation(m_samplesize, left, result, ::cos);
289        break;
290        case TAN:
291        tensor_unary_operation(m_samplesize, left, result, ::tan);
292        break;
293        case ASIN:
294        tensor_unary_operation(m_samplesize, left, result, ::asin);
295        break;
296        case ACOS:
297        tensor_unary_operation(m_samplesize, left, result, ::acos);
298        break;
299        case ATAN:
300        tensor_unary_operation(m_samplesize, left, result, ::atan);
301        break;
302        case SINH:
303        tensor_unary_operation(m_samplesize, left, result, ::sinh);
304        break;
305        case COSH:
306        tensor_unary_operation(m_samplesize, left, result, ::cosh);
307        break;
308        case TANH:
309        tensor_unary_operation(m_samplesize, left, result, ::tanh);
310        break;
311        case ERF:
312    #ifdef _WIN32
313        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
314    #else
315        tensor_unary_operation(m_samplesize, left, result, ::erf);
316        break;
317    #endif
318       case ASINH:
319    #ifdef _WIN32
320        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
321    #else
322        tensor_unary_operation(m_samplesize, left, result, ::asinh);
323    #endif  
324        break;
325       case ACOSH:
326    #ifdef _WIN32
327        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
328    #else
329        tensor_unary_operation(m_samplesize, left, result, ::acosh);
330    #endif  
331        break;
332       case ATANH:
333    #ifdef _WIN32
334        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
335    #else
336        tensor_unary_operation(m_samplesize, left, result, ::atanh);
337    #endif  
338        break;
339        case LOG10:
340        tensor_unary_operation(m_samplesize, left, result, ::log10);
341        break;
342        case LOG:
343        tensor_unary_operation(m_samplesize, left, result, ::log);
344        break;
345        case SIGN:
346        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
347        break;
348        case ABS:
349        tensor_unary_operation(m_samplesize, left, result, ::fabs);
350        break;
351        case NEG:
352        tensor_unary_operation(m_samplesize, left, result, negate<double>());
353        break;
354        case POS:
355        // it doesn't mean anything for delayed.
356        // it will just trigger a deep copy of the lazy object
357        throw DataException("Programmer error - POS not supported for lazy data.");
358        break;
359        case EXP:
360        tensor_unary_operation(m_samplesize, left, result, ::exp);
361        break;
362        case SQRT:
363        tensor_unary_operation(m_samplesize, left, result, ::sqrt);
364        break;
365        case RECIP:
366        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
367        break;
368        case GZ:
369        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
370        break;
371        case LZ:
372        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
373        break;
374        case GEZ:
375        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
376        break;
377        case LEZ:
378        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
379        break;
380    
381      default:      default:
382      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
383        }
384    }    }
385      return result;
386    }
387    
388    DataReady_ptr
389    DataLazy::resolve()
390    {
391      // This is broken!     We need to have a buffer per thread!
392      // so the allocation of v needs to move inside the loop somehow
393    
394    cout << "Sample size=" << m_samplesize << endl;
395    cout << "Buffers=" << m_buffsRequired << endl;
396    
397      size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);
398      int numthreads=1;
399    #ifdef _OPENMP
400      numthreads=omp_get_max_threads();
401      int threadnum=0;
402    #endif
403      ValueType v(numthreads*threadbuffersize);
404    cout << "Buffer created with size=" << v.size() << endl;
405      ValueType dummy(getNoValues());
406      DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);
407      ValueType& resvec=result->getVector();
408      DataReady_ptr resptr=DataReady_ptr(result);
409      int sample;
410      int resoffset;
411      int totalsamples=getNumSamples();
412      #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)
413      for (sample=0;sample<totalsamples;++sample)
414      {
415    #ifdef _OPENMP
416        const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());
417    #else
418        const double* res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.
419    #endif
420        resoffset=result->getPointOffset(sample,0);
421        for (int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector
422        {
423        resvec[resoffset]=res[i];
424        }
425      }
426      return resptr;
427  }  }
428    
429  std::string  std::string
430  DataLazy::toString() const  DataLazy::toString() const
431  {  {
432    return "Lazy evaluation object. No details available.";    ostringstream oss;
433      oss << "Lazy Data:";
434      intoString(oss);
435      return oss.str();
436    }
437    
438    void
439    DataLazy::intoString(ostringstream& oss) const
440    {
441      switch (getOpgroup(m_op))
442      {
443      case G_IDENTITY:
444        oss << '@' << m_id.get();
445        break;
446      case G_BINARY:
447        oss << '(';
448        m_left->intoString(oss);
449        oss << ' ' << opToString(m_op) << ' ';
450        m_right->intoString(oss);
451        oss << ')';
452        break;
453      case G_UNARY:
454        oss << opToString(m_op) << '(';
455        m_left->intoString(oss);
456        oss << ')';
457        break;
458      default:
459        oss << "UNKNOWN";
460      }
461  }  }
462    
463  // Note that in this case, deepCopy does not make copies of the leaves.  // Note that in this case, deepCopy does not make copies of the leaves.
# Line 160  DataLazy::deepCopy() Line 476  DataLazy::deepCopy()
476  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
477  DataLazy::getLength() const  DataLazy::getLength() const
478  {  {
479    return length;    return m_length;
480  }  }
481    
482    
483  DataAbstract*  DataAbstract*
484  DataLazy::getSlice(const DataTypes::RegionType& region) const  DataLazy::getSlice(const DataTypes::RegionType& region) const
485  {  {
486    // this seems like a really good one to include I just haven't added it yet    throw DataException("getSlice - not implemented for Lazy objects.");
   throw DataException("getSlice - not implemented for Lazy objects - yet.");  
487  }  }
488    
489  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type

Legend:
Removed from v.1868  
changed lines
  Added in v.1888

  ViewVC Help
Powered by ViewVC 1.1.26