/[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 854 by gross, Thu Sep 21 05:29:42 2006 UTC revision 2644 by jfenwick, Wed Sep 2 04:14:03 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  #ifndef PASO_MPI  #include "esysmpi.h"
 #define MPI_Comm long  
 #endif  
   
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 48  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 102  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 129  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 142  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);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \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);  
145    
146    /**    /**
147       \brief       \brief
148       Constructor which copies data from any object that can be converted into       Constructor which copies data from any object that can be treated like a python array/sequence.
      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 199  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 214  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 226  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       switches on update protection       switches on update protection
243    
244    */    */
245    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 248  class Data { Line 248  class Data {
248    
249    /**    /**
250       \brief       \brief
251       Returns trueif the data object is protected against update       Returns true, if the data object is protected against update
252    
253    */    */
254    ESCRIPT_DLL_API    ESCRIPT_DLL_API
255    bool    bool
256    isProtected() const;    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 values of all data-points 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    convertToNumArray();    setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
274    
275    /**    /**
276       \brief       \brief
277       Return the values of all data-points for the given sample as a single 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    const boost::python::numeric::array    void
281    convertToNumArrayFromSampleNo(int sampleNo);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282    
283    /**    /**
284       \brief       \brief
285       Return the value of the specified data-point as a single python numarray object.       sets the values of a data-point on this process
286    */    */
 #ifndef PASO_MPI    
   ESCRIPT_DLL_API  
   const boost::python::numeric::array  
   convertToNumArrayFromDPNo(int ProcNo,  
                                                         int sampleNo,  
                             int dataPointNo);  
 #else  
287    ESCRIPT_DLL_API    ESCRIPT_DLL_API
288    const boost::python::numeric::array    void
289    convertToNumArrayFromDPNo(int procNo,    setValueOfDataPoint(int dataPointNo, const double);
                 int sampleNo,  
                 int dataPointNo);  
 #endif  
   
290    
291    /**    /**
292       \brief       \brief Return a data point across all processors as a python tuple.
      Fills the expanded Data object from values of a python numarray object.  
293    */    */
294    ESCRIPT_DLL_API    ESCRIPT_DLL_API
295    void    const boost::python::object
296    fillFromNumArray(const boost::python::numeric::array);    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 316  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 325  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 363  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 379  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 394  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 427  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 447  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 477  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 507  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 568  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 593  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    
689    /**    /**
690       \brief       \brief
691       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
692       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
693       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
694         \param name - Input - name of tag.
695         \param value - Input - Value to associate with given key.
696      */
697      ESCRIPT_DLL_API
698      void
699      setTaggedValueByName(std::string name,
700                           const boost::python::object& value);
701    
702      /**
703         \brief
704         Assign the given value to the tag. Implicitly converts this
705         object to type DataTagged if it is constant.
706    
707       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
708       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
709      ==>*      ==>*
# Line 613  class Data { Line 716  class Data {
716    /**    /**
717       \brief       \brief
718       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
719       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
720       cannot be converted to a DataTagged object.  
721       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
722         \param pointshape - Input - The shape of the value parameter
723       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
724      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
725    */    */
726    ESCRIPT_DLL_API    ESCRIPT_DLL_API
727    void    void
728    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
729                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
730                            const DataTypes::ValueType& value,
731                int dataOffset=0);
732    
733    
734    
735    /**    /**
736      \brief      \brief
# Line 639  class Data { Line 747  class Data {
747    
748    /**    /**
749       \brief       \brief
750         set all values to zero
751         *
752      */
753      ESCRIPT_DLL_API
754      void
755      setToZero();
756    
757      /**
758         \brief
759       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
760       the result as a Data object.       the result as a Data object.
761       *       *
# Line 647  class Data { Line 764  class Data {
764    Data    Data
765    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
766    
767    
768      ESCRIPT_DLL_API
769      Data
770      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
771                           double undef, Data& B, double Bmin, double Bstep);
772    
773    
774      // This signature needs to be tuned
775      ESCRIPT_DLL_API
776      Data
777      interpolateFromTable(boost::python::object table, double Amin, double Astep,
778                           double undef, Data& B, double Bmin, double Bstep);
779    
780    /**    /**
781       \brief       \brief
782       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 662  class Data { Line 792  class Data {
792    grad() const;    grad() const;
793    
794    /**    /**
795       \brief      \brief
796       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
      *  
797    */    */
798    ESCRIPT_DLL_API    ESCRIPT_DLL_API
799    boost::python::numeric::array    boost::python::object
800    integrate() const;    integrateToTuple_const() const;
801    
802    
803      /**
804        \brief
805         Calculate the integral over the function space domain as a python tuple.
806      */
807      ESCRIPT_DLL_API
808      boost::python::object
809      integrateToTuple();
810    
811    
812    
813    /**    /**
814       \brief       \brief
# Line 735  class Data { Line 875  class Data {
875    /**    /**
876       \brief       \brief
877       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
878       *  
879         The method is not const because lazy data needs to be expanded before Lsup can be computed.
880         The _const form can be used when the Data object is const, however this will only work for
881         Data which is not Lazy.
882    
883         For Data which contain no samples (or tagged Data for which no tags in use have a value)
884         zero is returned.
885    */    */
886    ESCRIPT_DLL_API    ESCRIPT_DLL_API
887    double    double
888    Lsup() const;    Lsup();
889    
   /**  
      \brief  
      Return the minimum absolute value of this Data object.  
      *  
   */  
890    ESCRIPT_DLL_API    ESCRIPT_DLL_API
891    double    double
892    Linf() const;    Lsup_const() const;
893    
894    
895    /**    /**
896       \brief       \brief
897       Return the maximum value of this Data object.       Return the maximum value of this Data object.
898       *  
899         The method is not const because lazy data needs to be expanded before sup 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         a large negative value is returned.
905    */    */
906    ESCRIPT_DLL_API    ESCRIPT_DLL_API
907    double    double
908    sup() const;    sup();
909    
910      ESCRIPT_DLL_API
911      double
912      sup_const() const;
913    
914    
915    /**    /**
916       \brief       \brief
917       Return the minimum value of this Data object.       Return the minimum value of this Data object.
918       *  
919         The method is not const because lazy data needs to be expanded before inf 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 positive value is returned.
925    */    */
926    ESCRIPT_DLL_API    ESCRIPT_DLL_API
927    double    double
928    inf() const;    inf();
929    
930      ESCRIPT_DLL_API
931      double
932      inf_const() const;
933    
934    
935    
936    /**    /**
937       \brief       \brief
# Line 798  class Data { Line 963  class Data {
963    /**    /**
964       \brief       \brief
965       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
966       the minimum value in this Data object.       the minimum component value in this Data object.
967         \note If you are working in python, please consider using Locator
968    instead of manually manipulating process and point IDs.
969    */    */
970    ESCRIPT_DLL_API    ESCRIPT_DLL_API
971    const boost::python::tuple    const boost::python::tuple
972    mindp() const;    minGlobalDataPoint() const;
973    
974      /**
975         \brief
976         Return the (sample number, data-point number) of the data point with
977         the minimum component value in this Data object.
978         \note If you are working in python, please consider using Locator
979    instead of manually manipulating process and point IDs.
980      */
981    ESCRIPT_DLL_API    ESCRIPT_DLL_API
982    void    const boost::python::tuple
983    calc_mindp(int& ProcNo,    maxGlobalDataPoint() const;
984                          int& SampleNo,    
985               int& DataPointNo) const;  
986    
987    /**    /**
988       \brief       \brief
989       Return the sign of each data point of this Data object.       Return the sign of each data point of this Data object.
# Line 868  class Data { Line 1043  class Data {
1043    /**    /**
1044       \brief       \brief
1045       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.
1046       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
1047       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
1048       first non-zero entry is positive.       first non-zero entry is positive.
1049       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
1050       *       *
# Line 889  class Data { Line 1064  class Data {
1064    
1065    /**    /**
1066       \brief       \brief
1067         Return the error function erf of each data point of this Data object.
1068         *
1069      */
1070      ESCRIPT_DLL_API
1071      Data
1072      erf() const;
1073    
1074      /**
1075         \brief
1076       Return the sin of each data point of this Data object.       Return the sin of each data point of this Data object.
1077       *       *
1078    */    */
# Line 1064  class Data { Line 1248  class Data {
1248    /**    /**
1249       \brief       \brief
1250       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.
1251        
1252       \param right Input - the power to raise the object to.       \param right Input - the power to raise the object to.
1253       *       *
1254     */     */
# Line 1075  class Data { Line 1259  class Data {
1259    /**    /**
1260       \brief       \brief
1261       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.
1262        
1263       \param left Input - the bases       \param left Input - the bases
1264       *       *
1265     */     */
# Line 1100  class Data { Line 1284  class Data {
1284    void    void
1285    saveVTK(std::string fileName) const;    saveVTK(std::string fileName) const;
1286    
1287    
1288    
1289    /**    /**
1290       \brief       \brief
1291       Overloaded operator +=       Overloaded operator +=
# Line 1111  class Data { Line 1297  class Data {
1297    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1298    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1299    
1300      ESCRIPT_DLL_API
1301      Data& operator=(const Data& other);
1302    
1303    /**    /**
1304       \brief       \brief
1305       Overloaded operator -=       Overloaded operator -=
# Line 1203  class Data { Line 1392  class Data {
1392    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1393    inline    inline
1394    void    void
1395    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1396    
1397    /**    /**
1398       \brief       \brief
# Line 1214  class Data { Line 1403  class Data {
1403    */    */
1404    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1405    Data    Data
1406    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1407    
1408    /**    /**
1409       \brief       \brief
# Line 1227  class Data { Line 1416  class Data {
1416    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1417    void    void
1418    setSlice(const Data& value,    setSlice(const Data& value,
1419             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1420    
1421    /**    /**
1422       \brief       \brief
1423       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.  
1424    */    */
1425    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1426    void    void
1427    archiveData(const std::string fileName);          print(void);
1428    
1429    /**    /**
1430       \brief       \brief
1431       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1432       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1433       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.  
1434    */    */
1435    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1436    void          int
1437    extractData(const std::string fileName,          get_MPIRank(void) const;
               const FunctionSpace& fspace);  
   
1438    
1439    /**    /**
1440       \brief       \brief
1441       print the data values to stdout. Used for debugging       return the MPI rank number of the local data
1442                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1443                     is returned
1444    */    */
1445    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1446    void          int
1447      print(void);          get_MPISize(void) const;
1448    
1449    /**    /**
1450       \brief       \brief
1451       return the MPI rank number of the local data       return the MPI rank number of the local data
1452           MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()                   MPI_COMM_WORLD is assumed and returned.
          is returned  
1453    */    */
1454    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1455      int          MPI_Comm
1456      get_MPIRank(void) const;          get_MPIComm(void) const;
1457    
1458    /**    /**
1459       \brief       \brief
1460       return the MPI rank number of the local data       return the object produced by the factory, which is a DataConstant or DataExpanded
1461           MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()      TODO Ownership of this object should be explained in doco.
          is returned  
1462    */    */
1463    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1464      int          DataAbstract*
1465      get_MPISize(void) const;          borrowData(void) const;
1466    
1467      ESCRIPT_DLL_API
1468            DataAbstract_ptr
1469            borrowDataPtr(void) const;
1470    
   /**  
      \brief  
      return the MPI rank number of the local data  
          MPI_COMM_WORLD is assumed and returned.  
   */  
1471    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1472      MPI_Comm          DataReady_ptr
1473      get_MPIComm(void) const;          borrowReadyPtr(void) const;
1474    
1475    
1476    
1477    /**    /**
1478       \brief       \brief
1479       return the object produced by the factory, which is a DataConstant or DataExpanded       Return a pointer to the beginning of the datapoint at the specified offset.
1480         TODO Eventually these should be inlined.
1481         \param i - position(offset) in the underlying datastructure
1482    */    */
1483    
1484      ESCRIPT_DLL_API
1485            DataTypes::ValueType::const_reference
1486            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1487    
1488    
1489      ESCRIPT_DLL_API
1490            DataTypes::ValueType::reference
1491            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1492    
1493    
1494    
1495    /**
1496       \brief Create a buffer for use by getSample
1497       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1498       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1499      
1500       In multi-threaded sections, this needs to be called on each thread.
1501    
1502       \return A BufferGroup* if Data is lazy, NULL otherwise.
1503       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1504    */
1505    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1506      DataAbstract*    BufferGroup*
1507      borrowData(void) const;    allocSampleBuffer() const;
1508    
1509    /**
1510       \brief Free a buffer allocated with allocSampleBuffer.
1511       \param buffer Input - pointer to the buffer to deallocate.
1512    */
1513    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1514    
1515   protected:   protected:
1516    
1517   private:   private:
1518    
1519      double
1520      LsupWorker() const;
1521    
1522      double
1523      supWorker() const;
1524    
1525      double
1526      infWorker() const;
1527    
1528      boost::python::object
1529      integrateWorker() const;
1530    
1531      void
1532      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1533    
1534      void
1535      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1536    
1537    
1538    /**    /**
1539       \brief       \brief
1540       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1372  class Data { Line 1606  class Data {
1606       \brief       \brief
1607       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1608    */    */
1609    template <class IValueType>  
1610    void    void
1611    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1612             const DataTypes::ShapeType& shape,
1613               const FunctionSpace& what,               const FunctionSpace& what,
1614               bool expanded);               bool expanded);
1615    
1616      void
1617      initialise(const WrappedArray& value,
1618                     const FunctionSpace& what,
1619                     bool expanded);
1620    
1621    //    //
1622    // flag to protect the data object against any update    // flag to protect the data object against any update
1623    bool m_protected;    bool m_protected;
1624      mutable bool m_shared;
1625      bool m_lazy;
1626    
1627    //    //
1628    // pointer to the actual data object    // pointer to the actual data object
1629    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1630      DataAbstract_ptr m_data;
1631    
1632    // If possible please use getReadyPtr instead.
1633    // But see warning below.
1634      const DataReady*
1635      getReady() const;
1636    
1637      DataReady*
1638      getReady();
1639    
1640    
1641    // Be wary of using this for local operations since it (temporarily) increases reference count.
1642    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1643    // getReady() instead
1644      DataReady_ptr
1645      getReadyPtr();
1646    
1647      const_DataReady_ptr
1648      getReadyPtr() const;
1649    
1650    
1651      /**
1652       \brief Update the Data's shared flag
1653       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1654       For internal use only.
1655      */
1656      void updateShareStatus(bool nowshared) const
1657      {
1658        m_shared=nowshared;     // m_shared is mutable
1659      }
1660    
1661      // In the isShared() method below:
1662      // A problem would occur if m_data (the address pointed to) were being modified
1663      // while the call m_data->is_shared is being executed.
1664      //
1665      // Q: So why do I think this code can be thread safe/correct?
1666      // A: We need to make some assumptions.
1667      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1668      // 2. We assume that no constructions or assignments which will share previously unshared
1669      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1670    //    //
1671    // pointer to the internal profiling data    // This means that the only transition we need to consider, is when a previously shared object is
1672    struct profDataEntry *profData;    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1673      // In those cases the m_shared flag changes to false after m_data has completed changing.
1674      // For any threads executing before the flag switches they will assume the object is still shared.
1675      bool isShared() const
1676      {
1677        return m_shared;
1678    /*  if (m_shared) return true;
1679        if (m_data->isShared())        
1680        {                  
1681            updateShareStatus(true);
1682            return true;
1683        }
1684        return false;*/
1685      }
1686    
1687      void forceResolve()
1688      {
1689        if (isLazy())
1690        {
1691            #ifdef _OPENMP
1692            if (omp_in_parallel())
1693            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1694            throw DataException("Please do not call forceResolve() in a parallel region.");
1695            }
1696            #endif
1697            resolve();
1698        }
1699      }
1700    
1701      /**
1702      \brief if another object is sharing out member data make a copy to work with instead.
1703      This code should only be called from single threaded sections of code.
1704      */
1705      void exclusiveWrite()
1706      {
1707    #ifdef _OPENMP
1708      if (omp_in_parallel())
1709      {
1710    // *((int*)0)=17;
1711        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1712      }
1713    #endif
1714        forceResolve();
1715        if (isShared())
1716        {
1717            DataAbstract* t=m_data->deepCopy();
1718            set_m_data(DataAbstract_ptr(t));
1719        }
1720      }
1721    
1722      /**
1723      \brief checks if caller can have exclusive write to the object
1724      */
1725      void checkExclusiveWrite()
1726      {
1727        if  (isLazy() || isShared())
1728        {
1729            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1730        }
1731      }
1732    
1733      /**
1734      \brief Modify the data abstract hosted by this Data object
1735      For internal use only.
1736      Passing a pointer to null is permitted (do this in the destructor)
1737      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1738      */
1739      void set_m_data(DataAbstract_ptr p);
1740    
1741      friend class DataAbstract;        // To allow calls to updateShareStatus
1742    
1743  };  };
1744    
1745  template <class IValueType>  }   // end namespace escript
1746  void  
1747  Data::initialise(const IValueType& value,  
1748                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1749                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1750    // so that I can dynamic cast between them below.
1751    #include "DataReady.h"
1752    #include "DataLazy.h"
1753    
1754    namespace escript
1755  {  {
1756    //  
1757    // Construct a Data object of the appropriate type.  inline
1758    // Construct the object first as there seems to be a bug which causes  const DataReady*
1759    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1760    // within the shared_ptr constructor.  {
1761    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1762      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1763      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
1764      m_data=temp_data;  }
1765    } else {  
1766      DataAbstract* temp=new DataConstant(value,what);  inline
1767      boost::shared_ptr<DataAbstract> temp_data(temp);  DataReady*
1768      m_data=temp_data;  Data::getReady()
1769    }  {
1770       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1771       EsysAssert((dr!=0), "Error - casting to DataReady.");
1772       return dr;
1773    }
1774    
1775    // Be wary of using this for local operations since it (temporarily) increases reference count.
1776    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1777    // getReady() instead
1778    inline
1779    DataReady_ptr
1780    Data::getReadyPtr()
1781    {
1782       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1783       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1784       return dr;
1785  }  }
1786    
1787    
1788    inline
1789    const_DataReady_ptr
1790    Data::getReadyPtr() const
1791    {
1792       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1793       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1794       return dr;
1795    }
1796    
1797    inline
1798    DataAbstract::ValueType::value_type*
1799    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1800    {
1801       if (isLazy())
1802       {
1803        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1804       }
1805       return getReady()->getSampleDataRW(sampleNo);
1806    }
1807    
1808    inline
1809    const DataAbstract::ValueType::value_type*
1810    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1811    {
1812       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1813       if (l!=0)
1814       {
1815        size_t offset=0;
1816        if (bufferg==NULL)
1817        {
1818            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1819        }
1820        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1821        return &((*res)[offset]);
1822       }
1823       return getReady()->getSampleDataRO(sampleNo);
1824    }
1825    
1826    
1827    
1828    /**
1829       Modify a filename for MPI parallel output to multiple files
1830    */
1831    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1832    
1833  /**  /**
1834     Binary Data object operators.     Binary Data object operators.
1835  */  */
1836  inline double rpow(double x,double y)  inline double rpow(double x,double y)
1837  {  {
1838      return pow(y,x);      return pow(y,x);
1839  };  }
1840    
1841  /**  /**
1842    \brief    \brief
# Line 1514  ESCRIPT_DLL_API Data operator*(const boo Line 1930  ESCRIPT_DLL_API Data operator*(const boo
1930  */  */
1931  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);
1932    
1933    
1934    
1935  /**  /**
1936    \brief    \brief
1937    Output operator    Output operator
# Line 1523  ESCRIPT_DLL_API std::ostream& operator<< Line 1941  ESCRIPT_DLL_API std::ostream& operator<<
1941  /**  /**
1942    \brief    \brief
1943    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1944    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1945    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1946    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1947    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1948  */  */
1949  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1950  Data  Data
1951  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1952                       Data& arg1,                       Data& arg_1,
1953                       int axis_offset=0,                       int axis_offset=0,
1954                       int transpose=0);                       int transpose=0);
1955    
1956  /**  /**
1957    \brief    \brief
   Return true if operands are equivalent, else return false.  
   NB: this operator does very little at this point, and isn't to  
   be relied on. Requires further implementation.  
 */  
 //ESCRIPT_DLL_API bool operator==(const Data& left, const Data& right);  
   
 /**  
   \brief  
1958    Perform the given binary operation with this and right as operands.    Perform the given binary operation with this and right as operands.
1959    Right is a Data object.    Right is a Data object.
1960  */  */
# Line 1556  Data::binaryOp(const Data& right, Line 1966  Data::binaryOp(const Data& right,
1966  {  {
1967     //     //
1968     // 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
1969     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1970       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1971     }     }
1972    
1973       if (isLazy() || right.isLazy())
1974       {
1975         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1976       }
1977     //     //
1978     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1979     Data tempRight(right);     Data tempRight(right);
1980    
1981     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1982       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1983         //         //
1984         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1985         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1986       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1987         //         //
1988         // interpolate onto the RHS function space         // interpolate onto the RHS function space
1989         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
1990         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
1991           set_m_data(tempLeft.m_data);
1992       }       }
1993     }     }
1994     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1587  Data::binaryOp(const Data& right, Line 2004  Data::binaryOp(const Data& right,
2004       // of any data type       // of any data type
2005       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2006       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2007       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2008     } else if (isTagged()) {     } else if (isTagged()) {
2009       //       //
2010       // 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 1609  Data::binaryOp(const Data& right, Line 2026  Data::binaryOp(const Data& right,
2026       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2027       escript::binaryOp(*leftC,*rightC,operation);       escript::binaryOp(*leftC,*rightC,operation);
2028     }     }
    #if defined DOPROF  
    profData->binary++;  
    #endif  
 }  
   
 /**  
   \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);  
   }  
2029  }  }
2030    
2031  /**  /**
# Line 1684  Data::algorithm(BinaryFunction operation Line 2052  Data::algorithm(BinaryFunction operation
2052      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2053      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2054      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2055      } else if (isEmpty()) {
2056        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2057      } else if (isLazy()) {
2058        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2059      } else {
2060        throw DataException("Error - Data encapsulates an unknown type.");
2061    }    }
   return 0;  
2062  }  }
2063    
2064  /**  /**
2065    \brief    \brief
2066    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.
2067    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2068    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
2069    rank 0 Data object.    rank 0 Data object.
2070    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2071  */  */
# Line 1701  inline Line 2074  inline
2074  Data  Data
2075  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2076  {  {
2077    if (isExpanded()) {    if (isEmpty()) {
2078      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2079      }
2080      else if (isExpanded()) {
2081        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2082      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2083      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2084      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2085      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2086      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2087      return result;      return result;
2088    } else if (isTagged()) {    }
2089      else if (isTagged()) {
2090      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());  
2091      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2092      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2093        defval[0]=0;
2094        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2095      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2096      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2097    } else if (isConstant()) {    }
2098      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2099        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2100      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2101      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2102      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2103      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2104      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2105      return result;      return result;
2106      } else if (isLazy()) {
2107        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2108      } else {
2109        throw DataException("Error - Data encapsulates an unknown type.");
2110      }
2111    }
2112    
2113    /**
2114      \brief
2115      Compute a tensor operation with two Data objects
2116      \param arg_0 - Input - Data object
2117      \param arg_1 - Input - Data object
2118      \param operation - Input - Binary op functor
2119    */
2120    template <typename BinaryFunction>
2121    inline
2122    Data
2123    C_TensorBinaryOperation(Data const &arg_0,
2124                            Data const &arg_1,
2125                            BinaryFunction operation)
2126    {
2127      if (arg_0.isEmpty() || arg_1.isEmpty())
2128      {
2129         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2130      }
2131      if (arg_0.isLazy() || arg_1.isLazy())
2132      {
2133         throw DataException("Error - Operations not permitted on lazy data.");
2134      }
2135      // Interpolate if necessary and find an appropriate function space
2136      Data arg_0_Z, arg_1_Z;
2137      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2138        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2139          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2140          arg_1_Z = Data(arg_1);
2141        }
2142        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2143          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2144          arg_0_Z =Data(arg_0);
2145        }
2146        else {
2147          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2148        }
2149      } else {
2150          arg_0_Z = Data(arg_0);
2151          arg_1_Z = Data(arg_1);
2152      }
2153      // Get rank and shape of inputs
2154      int rank0 = arg_0_Z.getDataPointRank();
2155      int rank1 = arg_1_Z.getDataPointRank();
2156      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2157      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2158      int size0 = arg_0_Z.getDataPointSize();
2159      int size1 = arg_1_Z.getDataPointSize();
2160      // Declare output Data object
2161      Data res;
2162    
2163      if (shape0 == shape1) {
2164        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2165          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2166          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2167          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2168          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2169    
2170          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2171        }
2172        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2173    
2174          // Prepare the DataConstant input
2175          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2176    
2177          // Borrow DataTagged input from Data object
2178          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2179    
2180          // Prepare a DataTagged output 2
2181          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2182          res.tag();
2183          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2184    
2185          // Prepare offset into DataConstant
2186          int offset_0 = tmp_0->getPointOffset(0,0);
2187          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2188    
2189          // Get the pointers to the actual data
2190          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2191          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2192    
2193          // Compute a result for the default
2194          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2195          // Compute a result for each tag
2196          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2197          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2198          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2199            tmp_2->addTag(i->first);
2200            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2201            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2202    
2203            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2204          }
2205    
2206        }
2207        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2208          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2209          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2210          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2211          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2212    
2213          int sampleNo_1,dataPointNo_1;
2214          int numSamples_1 = arg_1_Z.getNumSamples();
2215          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2216          int offset_0 = tmp_0->getPointOffset(0,0);
2217          res.requireWrite();
2218          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2219          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2220            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2221              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2222              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2223              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2224              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2225              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2226              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2227            }
2228          }
2229    
2230        }
2231        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2232          // Borrow DataTagged input from Data object
2233          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2234    
2235          // Prepare the DataConstant input
2236          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2237    
2238          // Prepare a DataTagged output 2
2239          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2240          res.tag();
2241          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2242    
2243          // Prepare offset into DataConstant
2244          int offset_1 = tmp_1->getPointOffset(0,0);
2245    
2246          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2247          // Get the pointers to the actual data
2248          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2249          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2250          // Compute a result for the default
2251          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2252          // Compute a result for each tag
2253          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2254          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2255          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2256            tmp_2->addTag(i->first);
2257            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2258            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2259            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2260          }
2261    
2262        }
2263        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2264          // Borrow DataTagged input from Data object
2265          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2266    
2267          // Borrow DataTagged input from Data object
2268          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2269    
2270          // Prepare a DataTagged output 2
2271          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2272          res.tag();        // DataTagged output
2273          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2274    
2275          // Get the pointers to the actual data
2276          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2277          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2278          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2279    
2280          // Compute a result for the default
2281          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2282          // Merge the tags
2283          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2284          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2285          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2286          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2287            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2288          }
2289          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2290            tmp_2->addTag(i->first);
2291          }
2292          // Compute a result for each tag
2293          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2294          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2295    
2296            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2297            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2298            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2299    
2300            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2301          }
2302    
2303        }
2304        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2305          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2306          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2307          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2308          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2309          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2310    
2311          int sampleNo_0,dataPointNo_0;
2312          int numSamples_0 = arg_0_Z.getNumSamples();
2313          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2314          res.requireWrite();
2315          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2316          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2317            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2318            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2319            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2320              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2321              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2322              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2323              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2324              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2325            }
2326          }
2327    
2328        }
2329        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2330          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2331          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2332          DataConstant* tmp_1=dynamic_cast<DataConstant*>(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          int offset_1 = tmp_1->getPointOffset(0,0);
2339          res.requireWrite();
2340          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2341          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2342            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2343              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2344              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2345    
2346              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2347              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2348              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2349    
2350    
2351              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2352            }
2353          }
2354    
2355        }
2356        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2357          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2358          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2359          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2360          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2361          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2362    
2363          int sampleNo_0,dataPointNo_0;
2364          int numSamples_0 = arg_0_Z.getNumSamples();
2365          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2366          res.requireWrite();
2367          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2368          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2369            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2370            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2371            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2372              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2373              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2374              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2375              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2376              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2377            }
2378          }
2379    
2380        }
2381        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2382          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2383          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2384          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2385          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2386          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2387    
2388          int sampleNo_0,dataPointNo_0;
2389          int numSamples_0 = arg_0_Z.getNumSamples();
2390          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2391          res.requireWrite();
2392          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2393          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2394            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2395              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2396              int offset_1 = tmp_1->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              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2400              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2401              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2402            }
2403          }
2404    
2405        }
2406        else {
2407          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2408        }
2409    
2410      } else if (0 == rank0) {
2411        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2412          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2413          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2414          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2415          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2416          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2417        }
2418        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2419    
2420          // Prepare the DataConstant input
2421          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2422    
2423          // Borrow DataTagged input from Data object
2424          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2425    
2426          // Prepare a DataTagged output 2
2427          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2428          res.tag();
2429          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2430    
2431          // Prepare offset into DataConstant
2432          int offset_0 = tmp_0->getPointOffset(0,0);
2433          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2434    
2435          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2436          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2437    
2438          // Compute a result for the default
2439          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2440          // Compute a result for each tag
2441          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2442          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2443          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2444            tmp_2->addTag(i->first);
2445            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2446            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2447            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2448          }
2449    
2450        }
2451        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2452    
2453          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2454          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2455          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2456          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2457    
2458          int sampleNo_1,dataPointNo_1;
2459          int numSamples_1 = arg_1_Z.getNumSamples();
2460          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2461          int offset_0 = tmp_0->getPointOffset(0,0);
2462          res.requireWrite();
2463          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2464          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2465            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2466              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2467              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2468              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2469              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2470              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2471              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2472    
2473            }
2474          }
2475    
2476        }
2477        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2478    
2479          // Borrow DataTagged input from Data object
2480          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2481    
2482          // Prepare the DataConstant input
2483          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2484    
2485          // Prepare a DataTagged output 2
2486          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2487          res.tag();
2488          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2489    
2490          // Prepare offset into DataConstant
2491          int offset_1 = tmp_1->getPointOffset(0,0);
2492          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2493    
2494          // Get the pointers to the actual data
2495          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2496          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2497    
2498    
2499          // Compute a result for the default
2500          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2501          // Compute a result for each tag
2502          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2503          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2504          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2505            tmp_2->addTag(i->first);
2506            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2507            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2508    
2509            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2510          }
2511    
2512        }
2513        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2514    
2515          // Borrow DataTagged input from Data object
2516          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2517    
2518          // Borrow DataTagged input from Data object
2519          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2520    
2521          // Prepare a DataTagged output 2
2522          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2523          res.tag();        // DataTagged output
2524          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2525    
2526          // Get the pointers to the actual data
2527          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2528          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2529          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2530    
2531          // Compute a result for the default
2532          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2533          // Merge the tags
2534          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2535          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2536          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2537          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2538            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2539          }
2540          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2541            tmp_2->addTag(i->first);
2542          }
2543          // Compute a result for each tag
2544          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2545          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2546            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2547            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2548            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2549    
2550            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2551          }
2552    
2553        }
2554        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2555    
2556          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2557          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2558          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2559          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2560          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2561    
2562          int sampleNo_0,dataPointNo_0;
2563          int numSamples_0 = arg_0_Z.getNumSamples();
2564          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2565          res.requireWrite();
2566          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2567          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2568            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2569            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2570            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2571              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2572              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2573              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2574              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2575              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2576            }
2577          }
2578    
2579        }
2580        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2581          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2582          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2583          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2584          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2585    
2586          int sampleNo_0,dataPointNo_0;
2587          int numSamples_0 = arg_0_Z.getNumSamples();
2588          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2589          int offset_1 = tmp_1->getPointOffset(0,0);
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            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2594              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2595              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2596              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2597              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2598              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2599              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2600            }
2601          }
2602    
2603    
2604        }
2605        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2606    
2607          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2608          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2609          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2610          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2611          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2612    
2613          int sampleNo_0,dataPointNo_0;
2614          int numSamples_0 = arg_0_Z.getNumSamples();
2615          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2616          res.requireWrite();
2617          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2618          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2619            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2620            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2621            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2622              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2623              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2624              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2625              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2626              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2627            }
2628          }
2629    
2630        }
2631        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2632    
2633          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2634          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2635          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2636          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2637          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2638    
2639          int sampleNo_0,dataPointNo_0;
2640          int numSamples_0 = arg_0_Z.getNumSamples();
2641          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2642          res.requireWrite();
2643          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2644          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2645            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2646              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2647              int offset_1 = tmp_1->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              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2651              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2652              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2653            }
2654          }
2655    
2656        }
2657        else {
2658          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2659        }
2660    
2661      } else if (0 == rank1) {
2662        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2663          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2664          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2665          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2666          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2667          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2668        }
2669        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2670    
2671          // Prepare the DataConstant input
2672          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2673    
2674          // Borrow DataTagged input from Data object
2675          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2676    
2677          // Prepare a DataTagged output 2
2678          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2679          res.tag();
2680          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2681    
2682          // Prepare offset into DataConstant
2683          int offset_0 = tmp_0->getPointOffset(0,0);
2684          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2685    
2686          //Get the pointers to the actual data
2687          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2688          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2689    
2690          // Compute a result for the default
2691          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2692          // Compute a result for each tag
2693          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2694          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2695          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2696            tmp_2->addTag(i->first);
2697            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2698            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2699            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2700          }
2701        }
2702        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2703    
2704          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2705          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2706          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2707          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2708    
2709          int sampleNo_1,dataPointNo_1;
2710          int numSamples_1 = arg_1_Z.getNumSamples();
2711          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2712          int offset_0 = tmp_0->getPointOffset(0,0);
2713          res.requireWrite();
2714          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2715          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2716            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2717              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2718              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2719              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2720              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2721              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2722              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2723            }
2724          }
2725    
2726        }
2727        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2728    
2729          // Borrow DataTagged input from Data object
2730          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2731    
2732          // Prepare the DataConstant input
2733          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2734    
2735          // Prepare a DataTagged output 2
2736          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2737          res.tag();
2738          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2739    
2740          // Prepare offset into DataConstant
2741          int offset_1 = tmp_1->getPointOffset(0,0);
2742          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2743          // Get the pointers to the actual data
2744          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2745          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2746          // Compute a result for the default
2747          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2748          // Compute a result for each tag
2749          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2750          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2751          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2752            tmp_2->addTag(i->first);
2753            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2754            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2755            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2756          }
2757    
2758        }
2759        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2760    
2761          // Borrow DataTagged input from Data object
2762          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2763    
2764          // Borrow DataTagged input from Data object
2765          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2766    
2767          // Prepare a DataTagged output 2
2768          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2769          res.tag();        // DataTagged output
2770          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2771    
2772          // Get the pointers to the actual data
2773          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2774          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2775          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2776    
2777          // Compute a result for the default
2778          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2779          // Merge the tags
2780          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2781          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2782          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2783          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2784            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2785          }
2786          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2787            tmp_2->addTag(i->first);
2788          }
2789          // Compute a result for each tag
2790          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2791          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2792            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2793            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2794            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2795            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2796          }
2797    
2798        }
2799        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2800    
2801          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2802          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2803          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2804          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2805          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2806    
2807          int sampleNo_0,dataPointNo_0;
2808          int numSamples_0 = arg_0_Z.getNumSamples();
2809          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2810          res.requireWrite();
2811          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2812          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2813            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2814            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2815            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2816              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2817              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2818              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2819              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2820              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2821            }
2822          }
2823    
2824        }
2825        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2826          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2827          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2828          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2829          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2830    
2831          int sampleNo_0,dataPointNo_0;
2832          int numSamples_0 = arg_0_Z.getNumSamples();
2833          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2834          int offset_1 = tmp_1->getPointOffset(0,0);
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            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2839              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2840              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2841              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2842              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2843              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2844              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2845            }
2846          }
2847    
2848    
2849        }
2850        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2851    
2852          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2853          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2854          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2855          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2856          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2857    
2858          int sampleNo_0,dataPointNo_0;
2859          int numSamples_0 = arg_0_Z.getNumSamples();
2860          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2861          res.requireWrite();
2862          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2863          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2864            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2865            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2866            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2867              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2868              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2869              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2870              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2871              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2872            }
2873          }
2874    
2875        }
2876        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2877    
2878          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2879          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2880          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2881          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2882          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2883    
2884          int sampleNo_0,dataPointNo_0;
2885          int numSamples_0 = arg_0_Z.getNumSamples();
2886          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2887          res.requireWrite();
2888          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2889          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2890            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2891              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2892              int offset_1 = tmp_1->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              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2896              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2897              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2898            }
2899          }
2900    
2901        }
2902        else {
2903          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2904        }
2905    
2906      } else {
2907        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2908      }
2909    
2910      return res;
2911    }
2912    
2913    template <typename UnaryFunction>
2914    Data
2915    C_TensorUnaryOperation(Data const &arg_0,
2916                           UnaryFunction operation)
2917    {
2918      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2919      {
2920         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2921      }
2922      if (arg_0.isLazy())
2923      {
2924         throw DataException("Error - Operations not permitted on lazy data.");
2925      }
2926      // Interpolate if necessary and find an appropriate function space
2927      Data arg_0_Z = Data(arg_0);
2928    
2929      // Get rank and shape of inputs
2930      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2931      int size0 = arg_0_Z.getDataPointSize();
2932    
2933      // Declare output Data object
2934      Data res;
2935    
2936      if (arg_0_Z.isConstant()) {
2937        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2938        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2939        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2940        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2941      }
2942      else if (arg_0_Z.isTagged()) {
2943    
2944        // Borrow DataTagged input from Data object
2945        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2946    
2947        // Prepare a DataTagged output 2
2948        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2949        res.tag();
2950        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2951    
2952        // Get the pointers to the actual data
2953        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2954        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2955        // Compute a result for the default
2956        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2957        // Compute a result for each tag
2958        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2959        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2960        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2961          tmp_2->addTag(i->first);
2962          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2963          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2964          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2965        }
2966    
2967    }    }
2968    Data falseRetVal; // to keep compiler quiet    else if (arg_0_Z.isExpanded()) {
2969    return falseRetVal;  
2970        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2971        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2972        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2973    
2974        int sampleNo_0,dataPointNo_0;
2975        int numSamples_0 = arg_0_Z.getNumSamples();
2976        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2977        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2978        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2979          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2980            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2981            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2982            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2983            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2984            tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2985          }
2986        }
2987      }
2988      else {
2989        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
2990      }
2991    
2992      return res;
2993  }  }
2994    
2995  }  }

Legend:
Removed from v.854  
changed lines
  Added in v.2644

  ViewVC Help
Powered by ViewVC 1.1.26