/[escript]/branches/schroedinger/escript/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/schroedinger/escript/src/DataLazy.cpp

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

revision 1888 by jfenwick, Wed Oct 15 04:00:21 2008 UTC revision 1889 by jfenwick, Thu Oct 16 05:57:09 2008 UTC
# Line 123  calcBuffs(const DataLazy_ptr& left, cons Line 123  calcBuffs(const DataLazy_ptr& left, cons
123     }     }
124  }  }
125    
   
   
126  }   // end anonymous namespace  }   // end anonymous namespace
127    
128    
# Line 154  DataLazy::DataLazy(DataAbstract_ptr p) Line 152  DataLazy::DataLazy(DataAbstract_ptr p)
152     else     else
153     {     {
154      m_id=dynamic_pointer_cast<DataReady>(p);      m_id=dynamic_pointer_cast<DataReady>(p);
155        if(p->isConstant()) {m_readytype='C';}
156        else if(p->isExpanded()) {m_readytype='E';}
157        else if (p->isTagged()) {m_readytype='T';}
158        else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
159     }     }
160     m_length=p->getLength();     m_length=p->getLength();
161     m_buffsRequired=1;     m_buffsRequired=1;
# Line 178  DataLazy::DataLazy(DataAbstract_ptr left Line 180  DataLazy::DataLazy(DataAbstract_ptr left
180     {     {
181      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
182     }     }
183       m_readytype=lleft->m_readytype;
184     m_length=left->getLength();     m_length=left->getLength();
185     m_left=lleft;     m_left=lleft;
186     m_buffsRequired=1;     m_buffsRequired=1;
# Line 185  DataLazy::DataLazy(DataAbstract_ptr left Line 188  DataLazy::DataLazy(DataAbstract_ptr left
188  }  }
189    
190    
191  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
192      : parent(resultFS(left,right,op), resultShape(left,right,op)),  //  : parent(resultFS(left,right,op), resultShape(left,right,op)),
193      m_left(left),  //  m_left(left),
194      m_right(right),  //  m_right(right),
195      m_op(op)  //  m_op(op)
196  {  // {
197     if (getOpgroup(op)!=G_BINARY)  //    if (getOpgroup(op)!=G_BINARY)
198     {  //    {
199      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.");
200     }  //    }
201     m_length=resultLength(m_left,m_right,m_op);  //    m_length=resultLength(m_left,m_right,m_op);
202     m_samplesize=getNumDPPSample()*getNoValues();  //    m_samplesize=getNumDPPSample()*getNoValues();
203     m_buffsRequired=calcBuffs(m_left, m_right, m_op);  //    m_buffsRequired=calcBuffs(m_left, m_right, m_op);
204  cout << "(2)Lazy created with " << m_samplesize << endl;  // cout << "(2)Lazy created with " << m_samplesize << endl;
205  }  // }
206    
207  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
208      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
# Line 207  DataLazy::DataLazy(DataAbstract_ptr left Line 210  DataLazy::DataLazy(DataAbstract_ptr left
210  {  {
211     if (getOpgroup(op)!=G_BINARY)     if (getOpgroup(op)!=G_BINARY)
212     {     {
213      throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
214     }     }
215     if (left->isLazy())     if (left->isLazy())
216     {     {
# Line 225  DataLazy::DataLazy(DataAbstract_ptr left Line 228  DataLazy::DataLazy(DataAbstract_ptr left
228     {     {
229      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
230     }     }
231       char lt=m_left->m_readytype;
232       char rt=m_right->m_readytype;
233       if (lt=='E' || rt=='E')
234       {
235        m_readytype='E';
236       }
237       else if (lt=='T' || rt=='T')
238       {
239        m_readytype='T';
240       }
241       else
242       {
243        m_readytype='C';
244       }
245     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
246     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
247     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
# Line 245  DataLazy::getBuffsRequired() const Line 261  DataLazy::getBuffsRequired() const
261  }  }
262    
263    
264    DataReady_ptr
265    DataLazy::collapseToReady()
266    {
267      if (m_readytype=='E')
268      { // this is more an efficiency concern than anything else
269        throw DataException("Programmer Error - do not use collapse on Expanded data.");
270      }
271      if (m_op==IDENTITY)
272      {
273        return m_id;
274      }
275      DataReady_ptr pleft=m_left->collapseToReady();
276      Data left(pleft);
277      Data right;
278      if (getOpgroup(m_op)==G_BINARY)
279      {
280        right=Data(m_right->collapseToReady());
281      }
282      Data result;
283      switch(m_op)
284      {
285        case ADD:
286        result=left+right;
287        break;
288        case SUB:      
289        result=left-right;
290        break;
291        case MUL:      
292        result=left*right;
293        break;
294        case DIV:      
295        result=left/right;
296        break;
297        case SIN:
298        result=left.sin();  
299        break;
300        case COS:
301        result=left.cos();
302        break;
303        case TAN:
304        result=left.tan();
305        break;
306        case ASIN:
307        result=left.asin();
308        break;
309        case ACOS:
310        result=left.acos();
311        break;
312        case ATAN:
313        result=left.atan();
314        break;
315        case SINH:
316        result=left.sinh();
317        break;
318        case COSH:
319        result=left.cosh();
320        break;
321        case TANH:
322        result=left.tanh();
323        break;
324        case ERF:
325        result=left.erf();
326        break;
327       case ASINH:
328        result=left.asinh();
329        break;
330       case ACOSH:
331        result=left.acosh();
332        break;
333       case ATANH:
334        result=left.atanh();
335        break;
336        case LOG10:
337        result=left.log10();
338        break;
339        case LOG:
340        result=left.log();
341        break;
342        case SIGN:
343        result=left.sign();
344        break;
345        case ABS:
346        result=left.abs();
347        break;
348        case NEG:
349        result=left.neg();
350        break;
351        case POS:
352        // it doesn't mean anything for delayed.
353        // it will just trigger a deep copy of the lazy object
354        throw DataException("Programmer error - POS not supported for lazy data.");
355        break;
356        case EXP:
357        result=left.exp();
358        break;
359        case SQRT:
360        result=left.sqrt();
361        break;
362        case RECIP:
363        result=left.oneOver();
364        break;
365        case GZ:
366        result=left.wherePositive();
367        break;
368        case LZ:
369        result=left.whereNegative();
370        break;
371        case GEZ:
372        result=left.whereNonNegative();
373        break;
374        case LEZ:
375        result=left.whereNonPositive();
376        break;
377        default:
378        throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
379      }
380      return result.borrowReadyPtr();
381    }
382    
383    void
384    DataLazy::collapse()
385    {
386      if (m_op==IDENTITY)
387      {
388        return;
389      }
390      if (m_readytype=='E')
391      { // this is more an efficiency concern than anything else
392        throw DataException("Programmer Error - do not use collapse on Expanded data.");
393      }
394      m_id=collapseToReady();
395      m_op=IDENTITY;
396    }
397    
398    const double*
399    DataLazy::resolveUnary(ValueType& v,int sampleNo,  size_t offset) const
400    {
401        // we assume that any collapsing has been done before we get here
402        // since we only have one argument we don't need to think about only
403        // processing single points.
404      if (m_readytype!='E')
405      {
406        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
407      }
408      const double* left=m_left->resolveSample(v,sampleNo,offset);
409      double* result=&(v[offset]);
410      switch (m_op)
411      {
412        case SIN:  
413        tensor_unary_operation(m_samplesize, left, result, ::sin);
414        break;
415        case COS:
416        tensor_unary_operation(m_samplesize, left, result, ::cos);
417        break;
418        case TAN:
419        tensor_unary_operation(m_samplesize, left, result, ::tan);
420        break;
421        case ASIN:
422        tensor_unary_operation(m_samplesize, left, result, ::asin);
423        break;
424        case ACOS:
425        tensor_unary_operation(m_samplesize, left, result, ::acos);
426        break;
427        case ATAN:
428        tensor_unary_operation(m_samplesize, left, result, ::atan);
429        break;
430        case SINH:
431        tensor_unary_operation(m_samplesize, left, result, ::sinh);
432        break;
433        case COSH:
434        tensor_unary_operation(m_samplesize, left, result, ::cosh);
435        break;
436        case TANH:
437        tensor_unary_operation(m_samplesize, left, result, ::tanh);
438        break;
439        case ERF:
440    #ifdef _WIN32
441        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
442    #else
443        tensor_unary_operation(m_samplesize, left, result, ::erf);
444        break;
445    #endif
446       case ASINH:
447    #ifdef _WIN32
448        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
449    #else
450        tensor_unary_operation(m_samplesize, left, result, ::asinh);
451    #endif  
452        break;
453       case ACOSH:
454    #ifdef _WIN32
455        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
456    #else
457        tensor_unary_operation(m_samplesize, left, result, ::acosh);
458    #endif  
459        break;
460       case ATANH:
461    #ifdef _WIN32
462        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
463    #else
464        tensor_unary_operation(m_samplesize, left, result, ::atanh);
465    #endif  
466        break;
467        case LOG10:
468        tensor_unary_operation(m_samplesize, left, result, ::log10);
469        break;
470        case LOG:
471        tensor_unary_operation(m_samplesize, left, result, ::log);
472        break;
473        case SIGN:
474        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
475        break;
476        case ABS:
477        tensor_unary_operation(m_samplesize, left, result, ::fabs);
478        break;
479        case NEG:
480        tensor_unary_operation(m_samplesize, left, result, negate<double>());
481        break;
482        case POS:
483        // it doesn't mean anything for delayed.
484        // it will just trigger a deep copy of the lazy object
485        throw DataException("Programmer error - POS not supported for lazy data.");
486        break;
487        case EXP:
488        tensor_unary_operation(m_samplesize, left, result, ::exp);
489        break;
490        case SQRT:
491        tensor_unary_operation(m_samplesize, left, result, ::sqrt);
492        break;
493        case RECIP:
494        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
495        break;
496        case GZ:
497        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
498        break;
499        case LZ:
500        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
501        break;
502        case GEZ:
503        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
504        break;
505        case LEZ:
506        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
507        break;
508    
509        default:
510        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
511      }
512      return result;
513    }
514    
515    
516    #define PROC_OP(X) \
517        for (int i=0;i<steps;++i,resultp+=getNoValues()) \
518        { \
519    cout << "Step#" << i << " chunk=" << chunksize << endl; \
520    cout << left[0] << left[1] << left[2] << endl; \
521    cout << right[0] << right[1] << right[2] << endl; \
522           tensor_binary_operation(chunksize, left, right, resultp, X); \
523           left+=leftStep; \
524           right+=rightStep; \
525    cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \
526        }
527    
528    const double*
529    DataLazy::resolveBinary(ValueType& v,int sampleNo,  size_t offset) const
530    {
531        // again we assume that all collapsing has already been done
532        // so we have at least one expanded child.
533        // however, we could still have one of the children being not expanded.
534    
535    cout << "Resolve binary: " << toString() << endl;
536    
537      const double* left=m_left->resolveSample(v,sampleNo,offset);
538    cout << "Done Left " << left[0] << left[1] << left[2] << endl;
539      const double* right=m_right->resolveSample(v,sampleNo,offset);    
540    cout << "Done Right"  << right[0] << right[1] << right[2] << endl;
541        // now we need to know which args are expanded
542      bool leftExp=(m_left->m_readytype=='E');
543      bool rightExp=(m_right->m_readytype=='E');
544      bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step
545      int steps=(bigloops?1:getNumSamples());
546      size_t chunksize=(bigloops? m_samplesize : getNoValues());
547      int leftStep=(rightExp? getNoValues() : 0);
548      int rightStep=(leftExp? getNoValues() : 0);
549    cout << "left=" << left << " right=" << right << endl;
550      double* result=&(v[offset]);
551      double* resultp=result;
552      switch(m_op)
553      {
554        case ADD:
555        PROC_OP(plus<double>())
556        break;
557    // need to fill in the rest
558        default:
559        throw DataException("Programmer error - resolveBinay can not resolve operator "+opToString(m_op)+".");
560      }
561    cout << "About to return "  << result[0] << result[1] << result[2] << endl;;
562      return result;
563    }
564    
565    // the vector and the offset are a place where the method could write its data if it wishes
566    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
567    // Hence the return value to indicate where the data is actually stored.
568    // Regardless, the storage should be assumed to be used, even if it isn't.
569    const double*
570    DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset )
571    {
572    cout << "Resolve sample " << toString() << endl;
573        // collapse so we have a 'E' node or an IDENTITY for some other type
574      if (m_readytype!='E' && m_op!=IDENTITY)
575      {
576        collapse();
577      }
578    cout << "Post collapse check\n";
579      if (m_op==IDENTITY)  
580      {
581    cout << "In IDENTITY check\n";
582        const ValueType& vec=m_id->getVector();
583        if (m_readytype=='C')
584        {
585        return &(vec[0]);
586        }
587        return &(vec[m_id->getPointOffset(sampleNo, 0)]);
588      }
589      if (m_readytype!='E')
590      {
591        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
592      }
593    cout << "Calling sub resolvers\n";
594      switch (getOpgroup(m_op))
595      {
596      case G_UNARY: return resolveUnary(v,sampleNo,offset);
597      case G_BINARY: return resolveBinary(v,sampleNo,offset);
598      default:
599        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
600      }
601    }
602    
603    
604  // the vector and the offset are a place where the method could write its data if it wishes  // the vector and the offset are a place where the method could write its data if it wishes
605  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  // it is not obligated to do so. For example, if it has its own storage already, it can use that.
606  // Hence the return value to indicate where the data is actually stored.  // Hence the return value to indicate where the data is actually stored.
607  // Regardless, the storage should be assumed to be used, even if it isn't.  // Regardless, the storage should be assumed to be used, even if it isn't.
608  const double*  const double*
609  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  DataLazy::resolveSample2(ValueType& v,int sampleNo,  size_t offset )
610  {  {
611      if (m_readytype!='E')
612      {
613        throw DataException("Only supporting Expanded Data.");
614      }
615    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
616    {    {
617      const ValueType& vec=m_id->getVector();      const ValueType& vec=m_id->getVector();
# Line 388  DataLazy::resolveSample(ValueType& v,int Line 748  DataLazy::resolveSample(ValueType& v,int
748  DataReady_ptr  DataReady_ptr
749  DataLazy::resolve()  DataLazy::resolve()
750  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
751    
752  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
753  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
754    
755      if (m_readytype!='E')
756      {
757        collapse();
758      }
759      if (m_op==IDENTITY)
760      {
761        return m_id;
762      }
763        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
764    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);
765    int numthreads=1;    int numthreads=1;
766  #ifdef _OPENMP  #ifdef _OPENMP
767    numthreads=omp_get_max_threads();    numthreads=getNumberOfThreads();
768    int threadnum=0;    int threadnum=0;
769  #endif  #endif
770    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
# Line 418  cout << "Buffer created with size=" << v Line 785  cout << "Buffer created with size=" << v
785      const double* res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.      const double* res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.
786  #endif  #endif
787      resoffset=result->getPointOffset(sample,0);      resoffset=result->getPointOffset(sample,0);
788      for (int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector      for (unsigned int i=0;i<m_samplesize;++i,++resoffset)   // copy values into the output vector
789      {      {
790      resvec[resoffset]=res[i];      resvec[resoffset]=res[i];
791      }      }
# Line 440  DataLazy::intoString(ostringstream& oss) Line 807  DataLazy::intoString(ostringstream& oss)
807  {  {
808    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
809    {    {
810    case G_IDENTITY:    case G_IDENTITY:
811        if (m_id->isExpanded())
812        {
813           oss << "E";
814        }
815        else if (m_id->isTagged())
816        {
817          oss << "T";
818        }
819        else if (m_id->isConstant())
820        {
821          oss << "C";
822        }
823        else
824        {
825          oss << "?";
826        }
827      oss << '@' << m_id.get();      oss << '@' << m_id.get();
828      break;      break;
829    case G_BINARY:    case G_BINARY:

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

  ViewVC Help
Powered by ViewVC 1.1.26