/[escript]/trunk/escript/src/Data.h
ViewVC logotype

Diff of /trunk/escript/src/Data.h

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

revision 757 by woo409, Mon Jun 26 13:12:56 2006 UTC revision 2721 by jfenwick, Fri Oct 16 05:40:12 2009 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2009 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
8   *  Licensed under the Open Software License version 3.0    *  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Licensed under the Open Software License version 3.0
10   *                                                          *  * http://www.opensource.org/licenses/osl-3.0.php
11   ************************************************************  *
12  */  *******************************************************/
13    
14    
15  /** \file Data.h */  /** \file Data.h */
16    
# Line 17  Line 18 
18  #define DATA_H  #define DATA_H
19  #include "system_dep.h"  #include "system_dep.h"
20    
21    #include "DataTypes.h"
22  #include "DataAbstract.h"  #include "DataAbstract.h"
23  #include "DataAlgorithm.h"  #include "DataAlgorithm.h"
24  #include "FunctionSpace.h"  #include "FunctionSpace.h"
# Line 24  Line 26 
26  #include "UnaryOp.h"  #include "UnaryOp.h"
27  #include "DataException.h"  #include "DataException.h"
28    
29    
30  extern "C" {  extern "C" {
31  #include "DataC.h"  #include "DataC.h"
32  #include "paso/Paso.h"  //#include <omp.h>
33  }  }
34    
35    #include "esysmpi.h"
36  #include <string>  #include <string>
37  #include <algorithm>  #include <algorithm>
38    #include <sstream>
39    
40  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
41  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
42  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
43  #include <boost/python/numeric.hpp>  
44    #include "BufferGroup.h"
45    
46  namespace escript {  namespace escript {
47    
# Line 44  namespace escript { Line 50  namespace escript {
50  class DataConstant;  class DataConstant;
51  class DataTagged;  class DataTagged;
52  class DataExpanded;  class DataExpanded;
53    class DataLazy;
54    
55  /**  /**
56     \brief     \brief
57     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
58    
59     Description:     Description:
60     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
61     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
62     the object created for the lifetime of the object.     of the Data object.
63     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
64     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
65       Doing so will lead to invalid memory access.
66       This should not affect any methods exposed via boost::python.
67  */  */
68  class Data {  class Data {
69    
# Line 66  class Data { Line 74  class Data {
74    typedef double (*UnaryDFunPtr)(double);    typedef double (*UnaryDFunPtr)(double);
75    typedef double (*BinaryDFunPtr)(double,double);    typedef double (*BinaryDFunPtr)(double,double);
76    
77    
78    /**    /**
79       Constructors.       Constructors.
80    */    */
# Line 97  class Data { Line 106  class Data {
106         const FunctionSpace& what);         const FunctionSpace& what);
107    
108    /**    /**
109       \brief      \brief Copy Data from an existing vector
110       Constructor which copies data from a DataArrayView.    */
111    
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
112    ESCRIPT_DLL_API    ESCRIPT_DLL_API
113    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
114         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
115         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
116                     bool expanded=false);
117    
118    /**    /**
119       \brief       \brief
120       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
121    
122       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
123       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 124  class Data { Line 128  class Data {
128    */    */
129    ESCRIPT_DLL_API    ESCRIPT_DLL_API
130    Data(double value,    Data(double value,
131         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
132         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
133         bool expanded=false);         bool expanded=false);
134    
# Line 137  class Data { Line 141  class Data {
141    */    */
142    ESCRIPT_DLL_API    ESCRIPT_DLL_API
143    Data(const Data& inData,    Data(const Data& inData,
144         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Constructor which will create Tagged data if expanded is false.  
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
     ==>*  
   */  
   ESCRIPT_DLL_API  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
145    
146    /**    /**
147       \brief       \brief
148       Constructor which copies data from a python numarray.       Constructor which copies data from any object that can be treated like a python array/sequence.
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   ESCRIPT_DLL_API  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from any object that can be converted into  
      a python numarray.  
149    
150       \param value - Input - Input data.       \param value - Input - Input data.
151       \param what - Input - A description of what this data represents.       \param what - Input - A description of what this data represents.
# Line 194  class Data { Line 161  class Data {
161    /**    /**
162       \brief       \brief
163       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
164       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
165       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
166    
167       \param value - Input - Input data.       \param value - Input - Input data.
168       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
# Line 209  class Data { Line 176  class Data {
176       Constructor which creates a DataConstant of "shape" with constant value.       Constructor which creates a DataConstant of "shape" with constant value.
177    */    */
178    ESCRIPT_DLL_API    ESCRIPT_DLL_API
179    Data(double value,    Data(double value,
180         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
181         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
182         bool expanded=false);         bool expanded=false);
183    
184    
185    
186      /**
187        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
188        Once you have passed the pointer, do not delete it.
189      */
190      ESCRIPT_DLL_API
191      explicit Data(DataAbstract* underlyingdata);
192    
193      /**
194        \brief Create a Data based on the supplied DataAbstract
195      */
196      ESCRIPT_DLL_API
197      explicit Data(DataAbstract_ptr underlyingdata);
198    
199    /**    /**
200       \brief       \brief
201       Destructor       Destructor
# Line 221  class Data { Line 204  class Data {
204    ~Data();    ~Data();
205    
206    /**    /**
207       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
208    */    */
209    ESCRIPT_DLL_API    ESCRIPT_DLL_API
210    void    void
211    copy(const Data& other);    copy(const Data& other);
212    
213    /**    /**
214         \brief Return a pointer to a deep copy of this object.
215      */
216      ESCRIPT_DLL_API
217      Data
218      copySelf();
219    
220    
221      /**
222         \brief produce a delayed evaluation version of this Data.
223      */
224      ESCRIPT_DLL_API
225      Data
226      delay();
227    
228      /**
229         \brief convert the current data into lazy data.
230      */
231      ESCRIPT_DLL_API
232      void
233      delaySelf();
234    
235    
236      /**
237       Member access methods.       Member access methods.
238    */    */
239    
240    /**    /**
241       \brief       \brief
242       Return the values of all data-points as a single python numarray object.       switches on update protection
243    
244    */    */
245    ESCRIPT_DLL_API    ESCRIPT_DLL_API
246    const boost::python::numeric::array    void
247    convertToNumArray();    setProtection();
248    
249    /**    /**
250       \brief       \brief
251       Return the values of all data-points for the given sample as a single python numarray object.       Returns true, if the data object is protected against update
252    
253    */    */
254    ESCRIPT_DLL_API    ESCRIPT_DLL_API
255    const boost::python::numeric::array    bool
256    convertToNumArrayFromSampleNo(int sampleNo);    isProtected() const;
257    
258    
259    /**
260       \brief
261       Return the value of a data point as a python tuple.
262    */
263      ESCRIPT_DLL_API
264      const boost::python::object
265      getValueOfDataPointAsTuple(int dataPointNo);
266    
267    /**    /**
268       \brief       \brief
269       Return the value of the specified data-point as a single python numarray object.       sets the values of a data-point from a python object on this process
270    */    */
271    ESCRIPT_DLL_API    ESCRIPT_DLL_API
272    const boost::python::numeric::array    void
273    convertToNumArrayFromDPNo(int sampleNo,    setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
                             int dataPointNo);  
274    
275    /**    /**
276       \brief       \brief
277       Fills the expanded Data object from values of a python numarray object.       sets the values of a data-point from a array-like object on this process
278    */    */
279    ESCRIPT_DLL_API    ESCRIPT_DLL_API
280    void    void
281    fillFromNumArray(const boost::python::numeric::array);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282    
283      /**
284         \brief
285         sets the values of a data-point on this process
286      */
287      ESCRIPT_DLL_API
288      void
289      setValueOfDataPoint(int dataPointNo, const double);
290    
291      /**
292         \brief Return a data point across all processors as a python tuple.
293      */
294      ESCRIPT_DLL_API
295      const boost::python::object
296      getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
297    
298    /**    /**
299       \brief       \brief
300       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
301    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
302    */    */
303    ESCRIPT_DLL_API    ESCRIPT_DLL_API
304    int    int
# Line 284  class Data { Line 312  class Data {
312    escriptDataC    escriptDataC
313    getDataC();    getDataC();
314    
315    
316    
317    /**    /**
318       \brief       \brief
319       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 293  class Data { Line 323  class Data {
323    getDataC() const;    getDataC() const;
324    
325    /**    /**
326       \brief      \brief How much space is required to evaulate a sample of the Data.
      Write the data as a string.  
327    */    */
328    ESCRIPT_DLL_API    ESCRIPT_DLL_API
329    inline    size_t
330    std::string    getSampleBufferSize() const;
331    toString() const  
332    {  
     return m_data->toString();  
   }  
333    
334    /**    /**
335       \brief       \brief
336       Return the DataArrayView of the point data. This essentially contains       Write the data as a string. For large amounts of data, a summary is printed.
      the shape information for each data point although it also may be used  
      to manipulate the point data.  
337    */    */
338    ESCRIPT_DLL_API    ESCRIPT_DLL_API
339    inline    std::string
340    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
341    
342    /**    /**
343       \brief       \brief
# Line 331  class Data { Line 352  class Data {
352       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
353       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
354       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
355    */    */
356    ESCRIPT_DLL_API    ESCRIPT_DLL_API
357    void    void
358    tag();    tag();
359    
360    /**    /**
361        \brief If this data is lazy, then convert it to ready data.
362        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
363      */
364      ESCRIPT_DLL_API
365      void
366      resolve();
367    
368    
369      /**
370       \brief Ensures data is ready for write access.
371      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
372      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
373      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
374      Doing so might introduce additional sharing.
375      */
376      ESCRIPT_DLL_API
377      void
378      requireWrite();
379    
380      /**
381       \brief       \brief
382       Return true if this Data is expanded.       Return true if this Data is expanded.
383         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
384    */    */
385    ESCRIPT_DLL_API    ESCRIPT_DLL_API
386    bool    bool
# Line 347  class Data { Line 388  class Data {
388    
389    /**    /**
390       \brief       \brief
391         Return true if this Data is expanded or resolves to expanded.
392         That is, if it has a separate value for each datapoint in the sample.
393      */
394      ESCRIPT_DLL_API
395      bool
396      actsExpanded() const;
397      
398    
399      /**
400         \brief
401       Return true if this Data is tagged.       Return true if this Data is tagged.
402    */    */
403    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 362  class Data { Line 413  class Data {
413    isConstant() const;    isConstant() const;
414    
415    /**    /**
416         \brief Return true if this Data is lazy.
417      */
418      ESCRIPT_DLL_API
419      bool
420      isLazy() const;
421    
422      /**
423         \brief Return true if this data is ready.
424      */
425      ESCRIPT_DLL_API
426      bool
427      isReady() const;
428    
429      /**
430       \brief       \brief
431       Return true if this Data is empty.       Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
432    contains datapoints.
433    */    */
434    ESCRIPT_DLL_API    ESCRIPT_DLL_API
435    bool    bool
# Line 395  class Data { Line 461  class Data {
461    */    */
462    ESCRIPT_DLL_API    ESCRIPT_DLL_API
463    inline    inline
464    const AbstractDomain&  //   const AbstractDomain&
465      const_Domain_ptr
466    getDomain() const    getDomain() const
467    {    {
468       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
469    }    }
470    
471    
472      /**
473         \brief
474         Return the domain.
475         TODO: For internal use only.   This should be removed.
476      */
477      ESCRIPT_DLL_API
478      inline
479    //   const AbstractDomain&
480      Domain_ptr
481      getDomainPython() const
482      {
483         return getFunctionSpace().getDomainPython();
484      }
485    
486    /**    /**
487       \brief       \brief
488       Return a copy of the domain.       Return a copy of the domain.
# Line 415  class Data { Line 497  class Data {
497    */    */
498    ESCRIPT_DLL_API    ESCRIPT_DLL_API
499    inline    inline
500    int    unsigned int
501    getDataPointRank() const    getDataPointRank() const
502    {    {
503      return m_data->getPointDataView().getRank();      return m_data->getRank();
504    }    }
505    
506    /**    /**
507       \brief       \brief
508         Return the number of data points
509      */
510      ESCRIPT_DLL_API
511      inline
512      int
513      getNumDataPoints() const
514      {
515        return getNumSamples() * getNumDataPointsPerSample();
516      }
517      /**
518         \brief
519       Return the number of samples.       Return the number of samples.
520    */    */
521    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 445  class Data { Line 538  class Data {
538      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
539    }    }
540    
541    
542      /**
543        \brief
544        Return the number of values in the shape for this object.
545      */
546      ESCRIPT_DLL_API
547      int
548      getNoValues() const
549      {
550        return m_data->getNoValues();
551      }
552    
553    
554      /**
555         \brief
556         dumps the object into a netCDF file
557      */
558      ESCRIPT_DLL_API
559      void
560      dump(const std::string fileName) const;
561    
562     /**
563      \brief returns the values of the object as a list of tuples (one for each datapoint).
564    
565      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
566    If false, the result is a list of scalars [1, 2, ...]
567     */
568      ESCRIPT_DLL_API
569      const boost::python::object
570      toListOfTuples(bool scalarastuple=true);
571    
572    
573     /**
574        \brief
575        Return the sample data for the given sample no. This is not the
576        preferred interface but is provided for use by C code.
577        The bufferg parameter is only required for LazyData.
578        \param sampleNo - Input - the given sample no.
579        \param bufferg - A buffer to compute (and store) sample data in will be selected from this group.
580        \return pointer to the sample data.
581    */
582      ESCRIPT_DLL_API
583      inline
584      const DataAbstract::ValueType::value_type*
585      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg=0);
586    
587    
588    /**    /**
589       \brief       \brief
590       Return the sample data for the given sample no. This is not the       Return the sample data for the given sample no. This is not the
591       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
592       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
593         \return pointer to the sample data.
594    */    */
595    ESCRIPT_DLL_API    ESCRIPT_DLL_API
596    inline    inline
597    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
598    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
599    {  
     return m_data->getSampleData(sampleNo);  
   }  
600    
601    /**    /**
602       \brief       \brief
# Line 475  class Data { Line 614  class Data {
614    
615    /**    /**
616       \brief       \brief
617       Assign the given value to the data-points referenced by the given       Return a reference into the DataVector which points to the specified data point.
618       reference number.       \param sampleNo - Input -
619         \param dataPointNo - Input -
      The value supplied is a python numarray object.  The data from this numarray  
      is unpacked into a DataArray, and this is used to set the corresponding  
      data-points in the underlying Data object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Input - value to assign to data-points associated with  
                             the given reference number.  
620    */    */
621    ESCRIPT_DLL_API    ESCRIPT_DLL_API
622    void    DataTypes::ValueType::const_reference
623    setRefValue(int ref,    getDataPointRO(int sampleNo, int dataPointNo);
               const boost::python::numeric::array& value);  
624    
625    /**    /**
626       \brief       \brief
627       Return the values associated with the data-points referenced by the given       Return a reference into the DataVector which points to the specified data point.
628       reference number.       \param sampleNo - Input -
629         \param dataPointNo - Input -
      The value supplied is a python numarray object. The data from the corresponding  
      data-points in this Data object are packed into the given numarray object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Output - object to receive values from data-points  
                              associated with the given reference number.  
630    */    */
631    ESCRIPT_DLL_API    ESCRIPT_DLL_API
632    void    DataTypes::ValueType::reference
633    getRefValue(int ref,    getDataPointRW(int sampleNo, int dataPointNo);
634                boost::python::numeric::array& value);  
635    
636    
637    /**    /**
638       \brief       \brief
639       Return a view into the data for the data point specified.       Return the offset for the given sample and point within the sample
      NOTE: Construction of the DataArrayView is a relatively expensive  
      operation.  
      \param sampleNo - Input -  
      \param dataPointNo - Input -  
640    */    */
641    ESCRIPT_DLL_API    ESCRIPT_DLL_API
642    inline    inline
643    DataArrayView    DataTypes::ValueType::size_type
644    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
645                 int dataPointNo)                 int dataPointNo)
646    {    {
647      return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
648    }    }
649    
650    /**    /**
# Line 536  class Data { Line 652  class Data {
652       Return a reference to the data point shape.       Return a reference to the data point shape.
653    */    */
654    ESCRIPT_DLL_API    ESCRIPT_DLL_API
655    const DataArrayView::ShapeType&    inline
656    getDataPointShape() const;    const DataTypes::ShapeType&
657      getDataPointShape() const
658      {
659        return m_data->getShape();
660      }
661    
662    /**    /**
663       \brief       \brief
# Line 561  class Data { Line 681  class Data {
681       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
682    */    */
683    ESCRIPT_DLL_API    ESCRIPT_DLL_API
684    DataArrayView::ValueType::size_type    DataTypes::ValueType::size_type
685    getLength() const;    getLength() const;
686    
687    /**    /**
688      \brief Return true if this object contains no samples.
689      This is not the same as isEmpty()
690      */
691      ESCRIPT_DLL_API
692      bool
693      hasNoSamples() const
694      {
695        return getLength()==0;
696      }
697    
698      /**
699       \brief       \brief
700       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
701       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
702       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
703         \param name - Input - name of tag.
704         \param value - Input - Value to associate with given key.
705      */
706      ESCRIPT_DLL_API
707      void
708      setTaggedValueByName(std::string name,
709                           const boost::python::object& value);
710    
711      /**
712         \brief
713         Assign the given value to the tag. Implicitly converts this
714         object to type DataTagged if it is constant.
715    
716       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
717       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
718      ==>*      ==>*
# Line 581  class Data { Line 725  class Data {
725    /**    /**
726       \brief       \brief
727       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
728       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
729       cannot be converted to a DataTagged object.  
730       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
731         \param pointshape - Input - The shape of the value parameter
732       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
733      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
734    */    */
735    ESCRIPT_DLL_API    ESCRIPT_DLL_API
736    void    void
737    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
738                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
739                            const DataTypes::ValueType& value,
740                int dataOffset=0);
741    
742    
743    
744    /**    /**
745      \brief      \brief
# Line 607  class Data { Line 756  class Data {
756    
757    /**    /**
758       \brief       \brief
759         set all values to zero
760         *
761      */
762      ESCRIPT_DLL_API
763      void
764      setToZero();
765    
766      /**
767         \brief
768       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
769       the result as a Data object.       the result as a Data object.
770       *       *
# Line 615  class Data { Line 773  class Data {
773    Data    Data
774    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
775    
776    
777      ESCRIPT_DLL_API
778      Data
779      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
780                           double undef, Data& B, double Bmin, double Bstep,bool check_boundaries);
781    
782      ESCRIPT_DLL_API
783      Data
784      interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
785                           double undef,bool check_boundaries);
786    
787    
788    
789    
790      ESCRIPT_DLL_API
791      Data
792      interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
793                            Data& B, double Bmin, double Bstep, double undef,bool check_boundaries);
794    
795      ESCRIPT_DLL_API
796      Data
797      interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
798                            double undef,bool check_boundaries);
799    
800    /**    /**
801       \brief       \brief
802       Calculates the gradient of the data at the data points of functionspace.       Calculates the gradient of the data at the data points of functionspace.
# Line 630  class Data { Line 812  class Data {
812    grad() const;    grad() const;
813    
814    /**    /**
815        \brief
816         Calculate the integral over the function space domain as a python tuple.
817      */
818      ESCRIPT_DLL_API
819      boost::python::object
820      integrateToTuple_const() const;
821    
822    
823      /**
824        \brief
825         Calculate the integral over the function space domain as a python tuple.
826      */
827      ESCRIPT_DLL_API
828      boost::python::object
829      integrateToTuple();
830    
831    
832    
833      /**
834       \brief       \brief
835       Calculate the integral over the function space domain.       Returns 1./ Data object
836       *       *
837    */    */
838    ESCRIPT_DLL_API    ESCRIPT_DLL_API
839    boost::python::numeric::array    Data
840    integrate() const;    oneOver() const;
   
841    /**    /**
842       \brief       \brief
843       Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.       Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.
# Line 695  class Data { Line 895  class Data {
895    /**    /**
896       \brief       \brief
897       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
898       *  
899         The method is not const because lazy data needs to be expanded before Lsup can be computed.
900         The _const form can be used when the Data object is const, however this will only work for
901         Data which is not Lazy.
902    
903         For Data which contain no samples (or tagged Data for which no tags in use have a value)
904         zero is returned.
905    */    */
906    ESCRIPT_DLL_API    ESCRIPT_DLL_API
907    double    double
908    Lsup() const;    Lsup();
909    
   /**  
      \brief  
      Return the minimum absolute value of this Data object.  
      *  
   */  
910    ESCRIPT_DLL_API    ESCRIPT_DLL_API
911    double    double
912    Linf() const;    Lsup_const() const;
913    
914    
915    /**    /**
916       \brief       \brief
917       Return the maximum value of this Data object.       Return the maximum value of this Data object.
918       *  
919         The method is not const because lazy data needs to be expanded before sup can be computed.
920         The _const form can be used when the Data object is const, however this will only work for
921         Data which is not Lazy.
922    
923         For Data which contain no samples (or tagged Data for which no tags in use have a value)
924         a large negative value is returned.
925    */    */
926    ESCRIPT_DLL_API    ESCRIPT_DLL_API
927    double    double
928    sup() const;    sup();
929    
930      ESCRIPT_DLL_API
931      double
932      sup_const() const;
933    
934    
935    /**    /**
936       \brief       \brief
937       Return the minimum value of this Data object.       Return the minimum value of this Data object.
938       *  
939         The method is not const because lazy data needs to be expanded before inf can be computed.
940         The _const form can be used when the Data object is const, however this will only work for
941         Data which is not Lazy.
942    
943         For Data which contain no samples (or tagged Data for which no tags in use have a value)
944         a large positive value is returned.
945    */    */
946    ESCRIPT_DLL_API    ESCRIPT_DLL_API
947    double    double
948    inf() const;    inf();
949    
950      ESCRIPT_DLL_API
951      double
952      inf_const() const;
953    
954    
955    
956    /**    /**
957       \brief       \brief
# Line 758  class Data { Line 983  class Data {
983    /**    /**
984       \brief       \brief
985       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
986       the minimum value in this Data object.       the minimum component value in this Data object.
987         \note If you are working in python, please consider using Locator
988    instead of manually manipulating process and point IDs.
989    */    */
990    ESCRIPT_DLL_API    ESCRIPT_DLL_API
991    const boost::python::tuple    const boost::python::tuple
992    mindp() const;    minGlobalDataPoint() const;
993    
994      /**
995         \brief
996         Return the (sample number, data-point number) of the data point with
997         the minimum component value in this Data object.
998         \note If you are working in python, please consider using Locator
999    instead of manually manipulating process and point IDs.
1000      */
1001    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1002    void    const boost::python::tuple
1003    calc_mindp(int& SampleNo,    maxGlobalDataPoint() const;
1004               int& DataPointNo) const;  
1005    
1006    
1007    /**    /**
1008       \brief       \brief
# Line 781  class Data { Line 1016  class Data {
1016    
1017    /**    /**
1018       \brief       \brief
1019         Return the symmetric part of a matrix which is half the matrix plus its transpose.
1020         *
1021      */
1022      ESCRIPT_DLL_API
1023      Data
1024      symmetric() const;
1025    
1026      /**
1027         \brief
1028         Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
1029         *
1030      */
1031      ESCRIPT_DLL_API
1032      Data
1033      nonsymmetric() const;
1034    
1035      /**
1036         \brief
1037         Return the trace of a matrix
1038         *
1039      */
1040      ESCRIPT_DLL_API
1041      Data
1042      trace(int axis_offset) const;
1043    
1044      /**
1045         \brief
1046         Transpose each data point of this Data object around the given axis.
1047         *
1048      */
1049      ESCRIPT_DLL_API
1050      Data
1051      transpose(int axis_offset) const;
1052    
1053      /**
1054         \brief
1055       Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.       Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1056       Currently this function is restricted to rank 2, square shape, and dimension 3.       Currently this function is restricted to rank 2, square shape, and dimension 3.
1057       *       *
# Line 792  class Data { Line 1063  class Data {
1063    /**    /**
1064       \brief       \brief
1065       Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.       Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1066       the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than       the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1067       tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the       tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1068       first non-zero entry is positive.       first non-zero entry is positive.
1069       Currently this function is restricted to rank 2, square shape, and dimension 3       Currently this function is restricted to rank 2, square shape, and dimension 3
1070       *       *
# Line 804  class Data { Line 1075  class Data {
1075    
1076    /**    /**
1077       \brief       \brief
1078       Transpose each data point of this Data object around the given axis.       swaps the components axis0 and axis1
      --* not implemented yet *--  
1079       *       *
1080    */    */
1081    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1082    Data    Data
1083    transpose(int axis) const;    swapaxes(const int axis0, const int axis1) const;
1084    
1085    /**    /**
1086       \brief       \brief
1087       Calculate the trace of each data point of this Data object.       Return the error function erf of each data point of this Data object.
1088       *       *
1089    */    */
1090    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1091    Data    Data
1092    trace() const;    erf() const;
1093    
1094    /**    /**
1095       \brief       \brief
# Line 998  class Data { Line 1268  class Data {
1268    /**    /**
1269       \brief       \brief
1270       Return the given power of each data point of this boost python object.       Return the given power of each data point of this boost python object.
1271        
1272       \param right Input - the power to raise the object to.       \param right Input - the power to raise the object to.
1273       *       *
1274     */     */
# Line 1009  class Data { Line 1279  class Data {
1279    /**    /**
1280       \brief       \brief
1281       Return the given power of each data point of this boost python object.       Return the given power of each data point of this boost python object.
1282        
1283       \param left Input - the bases       \param left Input - the bases
1284       *       *
1285     */     */
# Line 1034  class Data { Line 1304  class Data {
1304    void    void
1305    saveVTK(std::string fileName) const;    saveVTK(std::string fileName) const;
1306    
1307    
1308    
1309    /**    /**
1310       \brief       \brief
1311       Overloaded operator +=       Overloaded operator +=
# Line 1045  class Data { Line 1317  class Data {
1317    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1318    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1319    
1320      ESCRIPT_DLL_API
1321      Data& operator=(const Data& other);
1322    
1323    /**    /**
1324       \brief       \brief
1325       Overloaded operator -=       Overloaded operator -=
# Line 1137  class Data { Line 1412  class Data {
1412    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1413    inline    inline
1414    void    void
1415    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1416    
1417    /**    /**
1418       \brief       \brief
# Line 1148  class Data { Line 1423  class Data {
1423    */    */
1424    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1425    Data    Data
1426    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1427    
1428    /**    /**
1429       \brief       \brief
# Line 1161  class Data { Line 1436  class Data {
1436    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1437    void    void
1438    setSlice(const Data& value,    setSlice(const Data& value,
1439             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1440    
1441    /**    /**
1442       \brief       \brief
1443       Archive the current Data object to the given file.       print the data values to stdout. Used for debugging
      \param fileName - Input - file to archive to.  
1444    */    */
1445    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1446    void    void
1447    archiveData(const std::string fileName);          print(void);
1448    
1449    /**    /**
1450       \brief       \brief
1451       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1452       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1453       Note - the current object must be of type DataEmpty.                   is returned
      \param fileName - Input - file to extract from.  
      \param fspace - Input - a suitable FunctionSpace descibing the data.  
1454    */    */
1455    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1456    void          int
1457    extractData(const std::string fileName,          get_MPIRank(void) const;
1458                const FunctionSpace& fspace);  
1459      /**
1460         \brief
1461         return the MPI rank number of the local data
1462                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1463                     is returned
1464      */
1465      ESCRIPT_DLL_API
1466            int
1467            get_MPISize(void) const;
1468    
1469      /**
1470         \brief
1471         return the MPI rank number of the local data
1472                     MPI_COMM_WORLD is assumed and returned.
1473      */
1474      ESCRIPT_DLL_API
1475            MPI_Comm
1476            get_MPIComm(void) const;
1477    
1478    /**    /**
1479       \brief       \brief
1480       print the data values to stdout. Used for debugging       return the object produced by the factory, which is a DataConstant or DataExpanded
1481        TODO Ownership of this object should be explained in doco.
1482    */    */
1483    void print();    ESCRIPT_DLL_API
1484            DataAbstract*
1485            borrowData(void) const;
1486    
1487      ESCRIPT_DLL_API
1488            DataAbstract_ptr
1489            borrowDataPtr(void) const;
1490    
1491      ESCRIPT_DLL_API
1492            DataReady_ptr
1493            borrowReadyPtr(void) const;
1494    
1495    
1496    
1497      /**
1498         \brief
1499         Return a pointer to the beginning of the datapoint at the specified offset.
1500         TODO Eventually these should be inlined.
1501         \param i - position(offset) in the underlying datastructure
1502      */
1503    
1504      ESCRIPT_DLL_API
1505            DataTypes::ValueType::const_reference
1506            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1507    
1508    
1509      ESCRIPT_DLL_API
1510            DataTypes::ValueType::reference
1511            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1512    
1513    
1514    
1515    /**
1516       \brief Create a buffer for use by getSample
1517       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1518       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1519      
1520       In multi-threaded sections, this needs to be called on each thread.
1521    
1522       \return A BufferGroup* if Data is lazy, NULL otherwise.
1523       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1524    */
1525      ESCRIPT_DLL_API
1526      BufferGroup*
1527      allocSampleBuffer() const;
1528    
1529    /**
1530       \brief Free a buffer allocated with allocSampleBuffer.
1531       \param buffer Input - pointer to the buffer to deallocate.
1532    */
1533    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1534    
1535   protected:   protected:
1536    
1537   private:   private:
1538    
1539    template <class BinaryOp>
1540      double
1541      lazyAlgWorker(double init, int mpiop_type);
1542    
1543      double
1544      LsupWorker() const;
1545    
1546      double
1547      supWorker() const;
1548    
1549      double
1550      infWorker() const;
1551    
1552      boost::python::object
1553      integrateWorker() const;
1554    
1555      void
1556      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1557    
1558      void
1559      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1560    
1561    
1562    /**    /**
1563       \brief       \brief
1564       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1249  class Data { Line 1612  class Data {
1612    
1613    /**    /**
1614       \brief       \brief
      Perform the given binary operation on all of the data's elements.  
      RHS is a boost::python object.  
   */  
   template <class BinaryFunction>  
   inline  
   void  
   binaryOp(const boost::python::object& right,  
            BinaryFunction operation);  
   
   /**  
      \brief  
1615       Convert the data type of the RHS to match this.       Convert the data type of the RHS to match this.
1616       \param right - Input - data type to match.       \param right - Input - data type to match.
1617    */    */
# Line 1278  class Data { Line 1630  class Data {
1630       \brief       \brief
1631       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1632    */    */
1633    template <class IValueType>  
1634    void    void
1635    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1636             const DataTypes::ShapeType& shape,
1637               const FunctionSpace& what,               const FunctionSpace& what,
1638               bool expanded);               bool expanded);
1639    
   /**  
      \brief  
      Reshape the data point if the data point is currently rank 0.  
      Will throw an exception if the data points are not rank 0.  
      The original data point value is used for all values of the new  
      data point.  
   */  
1640    void    void
1641    reshapeDataPoint(const DataArrayView::ShapeType& shape);    initialise(const WrappedArray& value,
1642                     const FunctionSpace& what,
1643                     bool expanded);
1644    
1645      //
1646      // flag to protect the data object against any update
1647      bool m_protected;
1648      mutable bool m_shared;
1649      bool m_lazy;
1650    
1651    //    //
1652    // pointer to the actual data object    // pointer to the actual data object
1653    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1654      DataAbstract_ptr m_data;
1655    
1656    // If possible please use getReadyPtr instead.
1657    // But see warning below.
1658      const DataReady*
1659      getReady() const;
1660    
1661      DataReady*
1662      getReady();
1663    
1664    
1665    // Be wary of using this for local operations since it (temporarily) increases reference count.
1666    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1667    // getReady() instead
1668      DataReady_ptr
1669      getReadyPtr();
1670    
1671      const_DataReady_ptr
1672      getReadyPtr() const;
1673    
1674    
1675      /**
1676       \brief Update the Data's shared flag
1677       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1678       For internal use only.
1679      */
1680      void updateShareStatus(bool nowshared) const
1681      {
1682        m_shared=nowshared;     // m_shared is mutable
1683      }
1684    
1685      // In the isShared() method below:
1686      // A problem would occur if m_data (the address pointed to) were being modified
1687      // while the call m_data->is_shared is being executed.
1688      //
1689      // Q: So why do I think this code can be thread safe/correct?
1690      // A: We need to make some assumptions.
1691      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1692      // 2. We assume that no constructions or assignments which will share previously unshared
1693      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1694    //    //
1695    // pointer to the internal profiling data    // This means that the only transition we need to consider, is when a previously shared object is
1696    struct profDataEntry *profData;    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1697      // In those cases the m_shared flag changes to false after m_data has completed changing.
1698      // For any threads executing before the flag switches they will assume the object is still shared.
1699      bool isShared() const
1700      {
1701        return m_shared;
1702    /*  if (m_shared) return true;
1703        if (m_data->isShared())        
1704        {                  
1705            updateShareStatus(true);
1706            return true;
1707        }
1708        return false;*/
1709      }
1710    
1711      void forceResolve()
1712      {
1713        if (isLazy())
1714        {
1715            #ifdef _OPENMP
1716            if (omp_in_parallel())
1717            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1718            throw DataException("Please do not call forceResolve() in a parallel region.");
1719            }
1720            #endif
1721            resolve();
1722        }
1723      }
1724    
1725      /**
1726      \brief if another object is sharing out member data make a copy to work with instead.
1727      This code should only be called from single threaded sections of code.
1728      */
1729      void exclusiveWrite()
1730      {
1731    #ifdef _OPENMP
1732      if (omp_in_parallel())
1733      {
1734    // *((int*)0)=17;
1735        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1736      }
1737    #endif
1738        forceResolve();
1739        if (isShared())
1740        {
1741            DataAbstract* t=m_data->deepCopy();
1742            set_m_data(DataAbstract_ptr(t));
1743        }
1744      }
1745    
1746      /**
1747      \brief checks if caller can have exclusive write to the object
1748      */
1749      void checkExclusiveWrite()
1750      {
1751        if  (isLazy() || isShared())
1752        {
1753            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1754        }
1755      }
1756    
1757      /**
1758      \brief Modify the data abstract hosted by this Data object
1759      For internal use only.
1760      Passing a pointer to null is permitted (do this in the destructor)
1761      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1762      */
1763      void set_m_data(DataAbstract_ptr p);
1764    
1765      friend class DataAbstract;        // To allow calls to updateShareStatus
1766    
1767  };  };
1768    
1769  template <class IValueType>  }   // end namespace escript
1770  void  
1771  Data::initialise(const IValueType& value,  
1772                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1773                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1774    // so that I can dynamic cast between them below.
1775    #include "DataReady.h"
1776    #include "DataLazy.h"
1777    
1778    namespace escript
1779  {  {
1780    //  
1781    // Construct a Data object of the appropriate type.  inline
1782    // Construct the object first as there seems to be a bug which causes  const DataReady*
1783    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1784    // within the shared_ptr constructor.  {
1785    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1786      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1787      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
1788      m_data=temp_data;  }
1789    } else {  
1790      DataAbstract* temp=new DataConstant(value,what);  inline
1791      boost::shared_ptr<DataAbstract> temp_data(temp);  DataReady*
1792      m_data=temp_data;  Data::getReady()
1793    }  {
1794       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1795       EsysAssert((dr!=0), "Error - casting to DataReady.");
1796       return dr;
1797    }
1798    
1799    // Be wary of using this for local operations since it (temporarily) increases reference count.
1800    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1801    // getReady() instead
1802    inline
1803    DataReady_ptr
1804    Data::getReadyPtr()
1805    {
1806       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1807       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1808       return dr;
1809    }
1810    
1811    
1812    inline
1813    const_DataReady_ptr
1814    Data::getReadyPtr() const
1815    {
1816       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1817       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1818       return dr;
1819    }
1820    
1821    inline
1822    DataAbstract::ValueType::value_type*
1823    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1824    {
1825       if (isLazy())
1826       {
1827        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1828       }
1829       return getReady()->getSampleDataRW(sampleNo);
1830    }
1831    
1832    inline
1833    const DataAbstract::ValueType::value_type*
1834    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1835    {
1836       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1837       if (l!=0)
1838       {
1839        size_t offset=0;
1840        if (bufferg==NULL)
1841        {
1842            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1843        }
1844        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1845        return &((*res)[offset]);
1846       }
1847       return getReady()->getSampleDataRO(sampleNo);
1848  }  }
1849    
1850    
1851    
1852    /**
1853       Modify a filename for MPI parallel output to multiple files
1854    */
1855    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1856    
1857  /**  /**
1858     Binary Data object operators.     Binary Data object operators.
1859  */  */
1860    inline double rpow(double x,double y)
1861    {
1862        return pow(y,x);
1863    }
1864    
1865  /**  /**
1866    \brief    \brief
# Line 1422  ESCRIPT_DLL_API Data operator*(const boo Line 1954  ESCRIPT_DLL_API Data operator*(const boo
1954  */  */
1955  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1956    
1957    
1958    
1959  /**  /**
1960    \brief    \brief
1961    Output operator    Output operator
# Line 1430  ESCRIPT_DLL_API std::ostream& operator<< Line 1964  ESCRIPT_DLL_API std::ostream& operator<<
1964    
1965  /**  /**
1966    \brief    \brief
1967    Return true if operands are equivalent, else return false.    Compute a tensor product of two Data objects
1968    NB: this operator does very little at this point, and isn't to    \param arg_0 - Input - Data object
1969    be relied on. Requires further implementation.    \param arg_1 - Input - Data object
1970      \param axis_offset - Input - axis offset
1971      \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1972  */  */
1973  //ESCRIPT_DLL_API bool operator==(const Data& left, const Data& right);  ESCRIPT_DLL_API
1974    Data
1975    C_GeneralTensorProduct(Data& arg_0,
1976                         Data& arg_1,
1977                         int axis_offset=0,
1978                         int transpose=0);
1979    
1980  /**  /**
1981    \brief    \brief
# Line 1449  Data::binaryOp(const Data& right, Line 1990  Data::binaryOp(const Data& right,
1990  {  {
1991     //     //
1992     // if this has a rank of zero promote it to the rank of the RHS     // if this has a rank of zero promote it to the rank of the RHS
1993     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1994       reshapeDataPoint(right.getPointDataView().getShape());       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1995       }
1996    
1997       if (isLazy() || right.isLazy())
1998       {
1999         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
2000     }     }
2001     //     //
2002     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
2003     Data tempRight(right);     Data tempRight(right);
2004    
2005     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
2006       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
2007         //         //
2008         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
2009         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
2010       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
2011         //         //
2012         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2013         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2014         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2015           set_m_data(tempLeft.m_data);
2016       }       }
2017     }     }
2018     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1480  Data::binaryOp(const Data& right, Line 2028  Data::binaryOp(const Data& right,
2028       // of any data type       // of any data type
2029       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2030       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2031       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2032     } else if (isTagged()) {     } else if (isTagged()) {
2033       //       //
2034       // Tagged data is operated on serially, the right hand side can be       // Tagged data is operated on serially, the right hand side can be
# Line 1506  Data::binaryOp(const Data& right, Line 2054  Data::binaryOp(const Data& right,
2054    
2055  /**  /**
2056    \brief    \brief
   Perform the given binary operation with this and right as operands.  
   Right is a boost::python object.  
 */  
 template <class BinaryFunction>  
 inline  
 void  
 Data::binaryOp(const boost::python::object& right,  
                BinaryFunction operation)  
 {  
    DataArray temp(right);  
    //  
    // if this has a rank of zero promote it to the rank of the RHS.  
    if (getPointDataView().getRank()==0 && temp.getView().getRank()!=0) {  
       reshapeDataPoint(temp.getView().getShape());  
    }  
    //  
    // Always allow scalar values for the RHS but check other shapes  
    if (temp.getView().getRank()!=0) {  
      if (!getPointDataView().checkShape(temp.getView().getShape())) {  
        throw DataException(getPointDataView().createShapeErrorMessage(  
                   "Error - RHS shape doesn't match LHS shape.",temp.getView().getShape()));  
      }  
    }  
    if (isExpanded()) {  
      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());  
      EsysAssert((leftC!=0),"Programming error - casting to DataExpanded.");  
      escript::binaryOp(*leftC,temp.getView(),operation);  
    } else if (isTagged()) {  
      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());  
      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");  
      escript::binaryOp(*leftC,temp.getView(),operation);  
    } else if (isConstant()) {  
      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());  
      EsysAssert((leftC!=0),"Programming error - casting to DataConstant.");  
      escript::binaryOp(*leftC,temp.getView(),operation);  
    }  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on other and return the result.  
   Given operation is performed on each element of each data point, thus  
   argument object is a rank n Data object, and returned object is a rank n  
   Data object.  
   Calls Data::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 Data  
 unaryOp(const Data& other,  
         UnaryFunction operation)  
 {  
   Data result;  
   result.copy(other);  
   result.unaryOp(operation);  
   return result;  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on this.  
   Given operation is performed on each element of each data point.  
   Calls escript::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 void  
 Data::unaryOp(UnaryFunction operation)  
 {  
   if (isExpanded()) {  
     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isTagged()) {  
     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isConstant()) {  
     DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");  
     escript::unaryOp(*leftC,operation);  
   }  
 }  
   
 /**  
   \brief  
2057    Perform the given Data object reduction algorithm on this and return the result.    Perform the given Data object reduction algorithm on this and return the result.
2058    Given operation combines each element of each data point, thus argument    Given operation combines each element of each data point, thus argument
2059    object (*this) is a rank n Data object, and returned object is a scalar.    object (*this) is a rank n Data object, and returned object is a scalar.
# Line 1614  Data::algorithm(BinaryFunction operation Line 2076  Data::algorithm(BinaryFunction operation
2076      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2077      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2078      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2079      } else if (isEmpty()) {
2080        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2081      } else if (isLazy()) {
2082        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2083      } else {
2084        throw DataException("Error - Data encapsulates an unknown type.");
2085    }    }
   return 0;  
2086  }  }
2087    
2088  /**  /**
2089    \brief    \brief
2090    Perform the given data point reduction algorithm on data and return the result.    Perform the given data point reduction algorithm on data and return the result.
2091    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2092    thus argument object is a rank n Data object, and returned object is a    thus argument object is a rank n Data object, and returned object is a
2093    rank 0 Data object.    rank 0 Data object.
2094    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2095  */  */
# Line 1631  inline Line 2098  inline
2098  Data  Data
2099  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2100  {  {
2101    if (isExpanded()) {    if (isEmpty()) {
2102      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2103      }
2104      else if (isExpanded()) {
2105        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2106      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2107      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2108      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2109      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2110      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2111      return result;      return result;
2112    } else if (isTagged()) {    }
2113      else if (isTagged()) {
2114      DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
     DataArrayView::ShapeType viewShape;  
     DataArrayView::ValueType viewData(1);  
     viewData[0]=0;  
     DataArrayView defaultValue(viewData,viewShape);  
     DataTagged::TagListType keys;  
     DataTagged::ValueListType values;  
     DataTagged::DataMapType::const_iterator i;  
     for (i=dataT->getTagLookup().begin();i!=dataT->getTagLookup().end();i++) {  
       keys.push_back(i->first);  
       values.push_back(defaultValue);  
     }  
     Data result(keys,values,defaultValue,getFunctionSpace());  
     DataTagged* resultT=dynamic_cast<DataTagged*>(result.m_data.get());  
2115      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2116      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2117        defval[0]=0;
2118        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2119      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2120      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2121    } else if (isConstant()) {    }
2122      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2123        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2124      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2125      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2126      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2127      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2128      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2129      return result;      return result;
2130      } else if (isLazy()) {
2131        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2132      } else {
2133        throw DataException("Error - Data encapsulates an unknown type.");
2134      }
2135    }
2136    
2137    /**
2138      \brief
2139      Compute a tensor operation with two Data objects
2140      \param arg_0 - Input - Data object
2141      \param arg_1 - Input - Data object
2142      \param operation - Input - Binary op functor
2143    */
2144    template <typename BinaryFunction>
2145    inline
2146    Data
2147    C_TensorBinaryOperation(Data const &arg_0,
2148                            Data const &arg_1,
2149                            BinaryFunction operation)
2150    {
2151      if (arg_0.isEmpty() || arg_1.isEmpty())
2152      {
2153         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2154      }
2155      if (arg_0.isLazy() || arg_1.isLazy())
2156      {
2157         throw DataException("Error - Operations not permitted on lazy data.");
2158      }
2159      // Interpolate if necessary and find an appropriate function space
2160      Data arg_0_Z, arg_1_Z;
2161      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2162        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2163          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2164          arg_1_Z = Data(arg_1);
2165        }
2166        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2167          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2168          arg_0_Z =Data(arg_0);
2169        }
2170        else {
2171          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2172        }
2173      } else {
2174          arg_0_Z = Data(arg_0);
2175          arg_1_Z = Data(arg_1);
2176      }
2177      // Get rank and shape of inputs
2178      int rank0 = arg_0_Z.getDataPointRank();
2179      int rank1 = arg_1_Z.getDataPointRank();
2180      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2181      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2182      int size0 = arg_0_Z.getDataPointSize();
2183      int size1 = arg_1_Z.getDataPointSize();
2184      // Declare output Data object
2185      Data res;
2186    
2187      if (shape0 == shape1) {
2188        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2189          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2190          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2191          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2192          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2193    
2194          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2195        }
2196        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2197    
2198          // Prepare the DataConstant input
2199          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2200    
2201          // Borrow DataTagged input from Data object
2202          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2203    
2204          // Prepare a DataTagged output 2
2205          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2206          res.tag();
2207          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2208    
2209          // Prepare offset into DataConstant
2210          int offset_0 = tmp_0->getPointOffset(0,0);
2211          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2212    
2213          // Get the pointers to the actual data
2214          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2215          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2216    
2217          // Compute a result for the default
2218          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2219          // Compute a result for each tag
2220          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2221          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2222          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2223            tmp_2->addTag(i->first);
2224            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2225            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2226    
2227            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2228          }
2229    
2230        }
2231        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2232          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2233          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2234          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2235          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2236    
2237          int sampleNo_1,dataPointNo_1;
2238          int numSamples_1 = arg_1_Z.getNumSamples();
2239          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2240          int offset_0 = tmp_0->getPointOffset(0,0);
2241          res.requireWrite();
2242          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2243          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2244            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2245              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2246              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2247              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2248              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2249              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2250              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2251            }
2252          }
2253    
2254        }
2255        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2256          // Borrow DataTagged input from Data object
2257          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2258    
2259          // Prepare the DataConstant input
2260          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2261    
2262          // Prepare a DataTagged output 2
2263          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2264          res.tag();
2265          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2266    
2267          // Prepare offset into DataConstant
2268          int offset_1 = tmp_1->getPointOffset(0,0);
2269    
2270          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2271          // Get the pointers to the actual data
2272          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2273          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2274          // Compute a result for the default
2275          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2276          // Compute a result for each tag
2277          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2278          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2279          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2280            tmp_2->addTag(i->first);
2281            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2282            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2283            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2284          }
2285    
2286        }
2287        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2288          // Borrow DataTagged input from Data object
2289          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2290    
2291          // Borrow DataTagged input from Data object
2292          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2293    
2294          // Prepare a DataTagged output 2
2295          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2296          res.tag();        // DataTagged output
2297          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2298    
2299          // Get the pointers to the actual data
2300          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2301          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2302          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2303    
2304          // Compute a result for the default
2305          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2306          // Merge the tags
2307          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2308          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2309          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2310          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2311            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2312          }
2313          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2314            tmp_2->addTag(i->first);
2315          }
2316          // Compute a result for each tag
2317          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2318          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2319    
2320            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2321            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2322            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2323    
2324            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2325          }
2326    
2327        }
2328        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2329          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2330          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2331          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2332          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2333          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2334    
2335          int sampleNo_0,dataPointNo_0;
2336          int numSamples_0 = arg_0_Z.getNumSamples();
2337          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2338          res.requireWrite();
2339          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2340          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2341            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2342            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2343            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2344              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2345              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2346              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2347              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2348              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2349            }
2350          }
2351    
2352        }
2353        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2354          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2355          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2356          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2357          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2358    
2359          int sampleNo_0,dataPointNo_0;
2360          int numSamples_0 = arg_0_Z.getNumSamples();
2361          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2362          int offset_1 = tmp_1->getPointOffset(0,0);
2363          res.requireWrite();
2364          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2365          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2366            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2367              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2368              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2369    
2370              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2371              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2372              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2373    
2374    
2375              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2376            }
2377          }
2378    
2379        }
2380        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2381          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2382          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2383          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2384          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2385          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2386    
2387          int sampleNo_0,dataPointNo_0;
2388          int numSamples_0 = arg_0_Z.getNumSamples();
2389          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2390          res.requireWrite();
2391          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2392          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2393            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2394            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2395            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2396              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2397              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2398              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2399              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2400              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2401            }
2402          }
2403    
2404        }
2405        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2406          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2407          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2408          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2409          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2410          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2411    
2412          int sampleNo_0,dataPointNo_0;
2413          int numSamples_0 = arg_0_Z.getNumSamples();
2414          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2415          res.requireWrite();
2416          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2417          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2418          dataPointNo_0=0;
2419    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2420              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2421              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2422              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2423              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2424              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2425              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2426              tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2427    //       }
2428          }
2429    
2430        }
2431        else {
2432          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2433        }
2434    
2435      } else if (0 == rank0) {
2436        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2437          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2438          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2439          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2440          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2441          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2442        }
2443        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2444    
2445          // Prepare the DataConstant input
2446          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2447    
2448          // Borrow DataTagged input from Data object
2449          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2450    
2451          // Prepare a DataTagged output 2
2452          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2453          res.tag();
2454          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2455    
2456          // Prepare offset into DataConstant
2457          int offset_0 = tmp_0->getPointOffset(0,0);
2458          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2459    
2460          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2461          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2462    
2463          // Compute a result for the default
2464          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2465          // Compute a result for each tag
2466          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2467          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2468          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2469            tmp_2->addTag(i->first);
2470            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2471            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2472            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2473          }
2474    
2475        }
2476        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2477    
2478          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2479          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2480          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2481          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2482    
2483          int sampleNo_1,dataPointNo_1;
2484          int numSamples_1 = arg_1_Z.getNumSamples();
2485          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2486          int offset_0 = tmp_0->getPointOffset(0,0);
2487          res.requireWrite();
2488          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2489          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2490            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2491              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2492              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2493              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2494              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2495              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2496              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2497    
2498            }
2499          }
2500    
2501        }
2502        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2503    
2504          // Borrow DataTagged input from Data object
2505          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2506    
2507          // Prepare the DataConstant input
2508          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2509    
2510          // Prepare a DataTagged output 2
2511          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2512          res.tag();
2513          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2514    
2515          // Prepare offset into DataConstant
2516          int offset_1 = tmp_1->getPointOffset(0,0);
2517          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2518    
2519          // Get the pointers to the actual data
2520          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2521          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2522    
2523    
2524          // Compute a result for the default
2525          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2526          // Compute a result for each tag
2527          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2528          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2529          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2530            tmp_2->addTag(i->first);
2531            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2532            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2533    
2534            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2535          }
2536    
2537        }
2538        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2539    
2540          // Borrow DataTagged input from Data object
2541          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2542    
2543          // Borrow DataTagged input from Data object
2544          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2545    
2546          // Prepare a DataTagged output 2
2547          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2548          res.tag();        // DataTagged output
2549          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2550    
2551          // Get the pointers to the actual data
2552          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2553          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2554          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2555    
2556          // Compute a result for the default
2557          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2558          // Merge the tags
2559          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2560          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2561          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2562          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2563            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2564          }
2565          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2566            tmp_2->addTag(i->first);
2567          }
2568          // Compute a result for each tag
2569          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2570          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2571            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2572            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2573            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2574    
2575            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2576          }
2577    
2578        }
2579        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2580    
2581          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2582          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2583          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2584          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2585          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2586    
2587          int sampleNo_0,dataPointNo_0;
2588          int numSamples_0 = arg_0_Z.getNumSamples();
2589          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2590          res.requireWrite();
2591          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2592          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2593            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2594            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2595            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2596              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2597              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2598              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2599              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2600              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2601            }
2602          }
2603    
2604        }
2605        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2606          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2607          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2608          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2609          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2610    
2611          int sampleNo_0,dataPointNo_0;
2612          int numSamples_0 = arg_0_Z.getNumSamples();
2613          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2614          int offset_1 = tmp_1->getPointOffset(0,0);
2615          res.requireWrite();
2616          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2617          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2618            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2619              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2620              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2621              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2622              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2623              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2624              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2625            }
2626          }
2627    
2628    
2629        }
2630        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2631    
2632          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2633          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2634          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2635          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2636          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2637    
2638          int sampleNo_0,dataPointNo_0;
2639          int numSamples_0 = arg_0_Z.getNumSamples();
2640          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2641          res.requireWrite();
2642          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2643          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2644            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2645            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2646            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2647              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2648              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2649              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2650              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2651              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2652            }
2653          }
2654    
2655        }
2656        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2657    
2658          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2659          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2660          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2661          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2662          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2663    
2664          int sampleNo_0,dataPointNo_0;
2665          int numSamples_0 = arg_0_Z.getNumSamples();
2666          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2667          res.requireWrite();
2668          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2669          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2670            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2671              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2672              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2673              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2674              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2675              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2676              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2677              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2678            }
2679          }
2680    
2681        }
2682        else {
2683          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2684        }
2685    
2686      } else if (0 == rank1) {
2687        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2688          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2689          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2690          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2691          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2692          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2693        }
2694        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2695    
2696          // Prepare the DataConstant input
2697          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2698    
2699          // Borrow DataTagged input from Data object
2700          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2701    
2702          // Prepare a DataTagged output 2
2703          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2704          res.tag();
2705          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2706    
2707          // Prepare offset into DataConstant
2708          int offset_0 = tmp_0->getPointOffset(0,0);
2709          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2710    
2711          //Get the pointers to the actual data
2712          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2713          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2714    
2715          // Compute a result for the default
2716          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2717          // Compute a result for each tag
2718          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2719          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2720          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2721            tmp_2->addTag(i->first);
2722            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2723            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2724            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2725          }
2726        }
2727        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2728    
2729          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2730          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2731          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2732          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2733    
2734          int sampleNo_1,dataPointNo_1;
2735          int numSamples_1 = arg_1_Z.getNumSamples();
2736          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2737          int offset_0 = tmp_0->getPointOffset(0,0);
2738          res.requireWrite();
2739          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2740          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2741            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2742              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2743              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2744              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2745              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2746              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2747              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2748            }
2749          }
2750    
2751        }
2752        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2753    
2754          // Borrow DataTagged input from Data object
2755          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2756    
2757          // Prepare the DataConstant input
2758          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2759    
2760          // Prepare a DataTagged output 2
2761          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2762          res.tag();
2763          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2764    
2765          // Prepare offset into DataConstant
2766          int offset_1 = tmp_1->getPointOffset(0,0);
2767          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2768          // Get the pointers to the actual data
2769          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2770          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2771          // Compute a result for the default
2772          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2773          // Compute a result for each tag
2774          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2775          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2776          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2777            tmp_2->addTag(i->first);
2778            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2779            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2780            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2781          }
2782    
2783        }
2784        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2785    
2786          // Borrow DataTagged input from Data object
2787          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2788    
2789          // Borrow DataTagged input from Data object
2790          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2791    
2792          // Prepare a DataTagged output 2
2793          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2794          res.tag();        // DataTagged output
2795          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2796    
2797          // Get the pointers to the actual data
2798          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2799          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2800          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2801    
2802          // Compute a result for the default
2803          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2804          // Merge the tags
2805          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2806          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2807          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2808          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2809            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2810          }
2811          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2812            tmp_2->addTag(i->first);
2813          }
2814          // Compute a result for each tag
2815          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2816          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2817            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2818            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2819            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2820            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2821          }
2822    
2823        }
2824        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2825    
2826          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2827          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2828          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2829          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2830          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2831    
2832          int sampleNo_0,dataPointNo_0;
2833          int numSamples_0 = arg_0_Z.getNumSamples();
2834          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2835          res.requireWrite();
2836          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2837          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2838            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2839            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2840            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2841              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2842              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2843              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2844              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2845              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2846            }
2847          }
2848    
2849        }
2850        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2851          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2852          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2853          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2854          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2855    
2856          int sampleNo_0,dataPointNo_0;
2857          int numSamples_0 = arg_0_Z.getNumSamples();
2858          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2859          int offset_1 = tmp_1->getPointOffset(0,0);
2860          res.requireWrite();
2861          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2862          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2863            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2864              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2865              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2866              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2867              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2868              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2869              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2870            }
2871          }
2872    
2873    
2874        }
2875        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2876    
2877          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2878          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2879          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2880          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2881          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2882    
2883          int sampleNo_0,dataPointNo_0;
2884          int numSamples_0 = arg_0_Z.getNumSamples();
2885          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2886          res.requireWrite();
2887          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2888          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2889            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2890            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2891            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2892              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2893              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2894              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2895              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2896              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2897            }
2898          }
2899    
2900        }
2901        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2902    
2903          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2904          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2905          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2906          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2907          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2908    
2909          int sampleNo_0,dataPointNo_0;
2910          int numSamples_0 = arg_0_Z.getNumSamples();
2911          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2912          res.requireWrite();
2913          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2914          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2915            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2916              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2917              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2918              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2919              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2920              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2921              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2922              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2923            }
2924          }
2925    
2926        }
2927        else {
2928          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2929        }
2930    
2931      } else {
2932        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2933    }    }
2934    Data falseRetVal; // to keep compiler quiet  
2935    return falseRetVal;    return res;
2936    }
2937    
2938    template <typename UnaryFunction>
2939    Data
2940    C_TensorUnaryOperation(Data const &arg_0,
2941                           UnaryFunction operation)
2942    {
2943      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2944      {
2945         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2946      }
2947      if (arg_0.isLazy())
2948      {
2949         throw DataException("Error - Operations not permitted on lazy data.");
2950      }
2951      // Interpolate if necessary and find an appropriate function space
2952      Data arg_0_Z = Data(arg_0);
2953    
2954      // Get rank and shape of inputs
2955      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2956      int size0 = arg_0_Z.getDataPointSize();
2957    
2958      // Declare output Data object
2959      Data res;
2960    
2961      if (arg_0_Z.isConstant()) {
2962        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2963        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2964        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2965        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2966      }
2967      else if (arg_0_Z.isTagged()) {
2968    
2969        // Borrow DataTagged input from Data object
2970        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2971    
2972        // Prepare a DataTagged output 2
2973        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2974        res.tag();
2975        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2976    
2977        // Get the pointers to the actual data
2978        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2979        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2980        // Compute a result for the default
2981        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2982        // Compute a result for each tag
2983        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2984        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2985        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2986          tmp_2->addTag(i->first);
2987          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2988          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2989          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2990        }
2991    
2992      }
2993      else if (arg_0_Z.isExpanded()) {
2994    
2995        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2996        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2997        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2998    
2999        int sampleNo_0,dataPointNo_0;
3000        int numSamples_0 = arg_0_Z.getNumSamples();
3001        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3002        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3003        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3004        dataPointNo_0=0;
3005    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3006            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3007            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3008            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3009            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3010            tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3011    //      }
3012        }
3013      }
3014      else {
3015        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3016      }
3017    
3018      return res;
3019  }  }
3020    
3021  }  }

Legend:
Removed from v.757  
changed lines
  Added in v.2721

  ViewVC Help
Powered by ViewVC 1.1.26