/[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 922 by gross, Fri Jan 5 04:23:05 2007 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 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;
   /**  
      \brief  
      Return the values of all data-points as a single python numarray object.  
   */  
   ESCRIPT_DLL_API  
   const boost::python::numeric::array  
   convertToNumArray();  
257    
258    /**  
259       \brief  /**
260       Fills the expanded Data object from values of a python numarray object.     \brief
261    */     Return the value of a data point as a python tuple.
262    */
263    ESCRIPT_DLL_API    ESCRIPT_DLL_API
264    void    const boost::python::object
265    fillFromNumArray(const boost::python::numeric::array);    getValueOfDataPointAsTuple(int dataPointNo);
266    
267    /**    /**
268       \brief       \brief
269       Return the values of a data point on this process       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    getValueOfDataPoint(int dataPointNo);    setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
274    
275    /**    /**
276       \brief       \brief
277       sets the values of a data-point on this process       sets the values of a data-point from a array-like object on this process
278    */    */
279    ESCRIPT_DLL_API    ESCRIPT_DLL_API
280    void    void
281    setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282    
283    /**    /**
284       \brief       \brief
# Line 295  class Data { Line 289  class Data {
289    setValueOfDataPoint(int dataPointNo, const double);    setValueOfDataPoint(int dataPointNo, const double);
290    
291    /**    /**
292       \brief       \brief Return a data point across all processors as a python tuple.
      Return the value of the specified data-point across all processors  
293    */    */
294    ESCRIPT_DLL_API    ESCRIPT_DLL_API
295    const boost::python::numeric::array    const boost::python::object
296    getValueOfGlobalDataPoint(int procNo, int dataPointNo);    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 321  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 330  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 368  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 384  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 399  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 432  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 452  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    /**    /**
# Line 493  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 523  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 584  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 609  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 629  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 655  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 662  class Data { Line 763  class Data {
763    ESCRIPT_DLL_API    ESCRIPT_DLL_API
764    Data    Data
765    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
   
766    /**    /**
767       \brief       \brief
768       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 678  class Data { Line 778  class Data {
778    grad() const;    grad() const;
779    
780    /**    /**
781       \brief      \brief
782       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
783       *    */
784      ESCRIPT_DLL_API
785      boost::python::object
786      integrateToTuple_const() const;
787    
788    
789      /**
790        \brief
791         Calculate the integral over the function space domain as a python tuple.
792    */    */
793    ESCRIPT_DLL_API    ESCRIPT_DLL_API
794    boost::python::numeric::array    boost::python::object
795    integrate() const;    integrateToTuple();
796    
797    
798    
799    /**    /**
800       \brief       \brief
# Line 751  class Data { Line 861  class Data {
861    /**    /**
862       \brief       \brief
863       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
864       *  
865         The method is not const because lazy data needs to be expanded before Lsup can be computed.
866         The _const form can be used when the Data object is const, however this will only work for
867         Data which is not Lazy.
868    
869         For Data which contain no samples (or tagged Data for which no tags in use have a value)
870         zero is returned.
871    */    */
872    ESCRIPT_DLL_API    ESCRIPT_DLL_API
873    double    double
874    Lsup() const;    Lsup();
875    
   /**  
      \brief  
      Return the minimum absolute value of this Data object.  
      *  
   */  
876    ESCRIPT_DLL_API    ESCRIPT_DLL_API
877    double    double
878    Linf() const;    Lsup_const() const;
879    
880    
881    /**    /**
882       \brief       \brief
883       Return the maximum value of this Data object.       Return the maximum value of this Data object.
884       *  
885         The method is not const because lazy data needs to be expanded before sup can be computed.
886         The _const form can be used when the Data object is const, however this will only work for
887         Data which is not Lazy.
888    
889         For Data which contain no samples (or tagged Data for which no tags in use have a value)
890         a large negative value is returned.
891    */    */
892    ESCRIPT_DLL_API    ESCRIPT_DLL_API
893    double    double
894    sup() const;    sup();
895    
896      ESCRIPT_DLL_API
897      double
898      sup_const() const;
899    
900    
901    /**    /**
902       \brief       \brief
903       Return the minimum value of this Data object.       Return the minimum value of this Data object.
904       *  
905         The method is not const because lazy data needs to be expanded before inf can be computed.
906         The _const form can be used when the Data object is const, however this will only work for
907         Data which is not Lazy.
908    
909         For Data which contain no samples (or tagged Data for which no tags in use have a value)
910         a large positive value is returned.
911    */    */
912    ESCRIPT_DLL_API    ESCRIPT_DLL_API
913    double    double
914    inf() const;    inf();
915    
916      ESCRIPT_DLL_API
917      double
918      inf_const() const;
919    
920    
921    
922    /**    /**
923       \brief       \brief
# Line 814  class Data { Line 949  class Data {
949    /**    /**
950       \brief       \brief
951       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
952       the minimum value in this Data object.       the minimum component value in this Data object.
953         \note If you are working in python, please consider using Locator
954    instead of manually manipulating process and point IDs.
955    */    */
956    ESCRIPT_DLL_API    ESCRIPT_DLL_API
957    const boost::python::tuple    const boost::python::tuple
958    minGlobalDataPoint() const;    minGlobalDataPoint() const;
959    
960      /**
961         \brief
962         Return the (sample number, data-point number) of the data point with
963         the minimum component value in this Data object.
964         \note If you are working in python, please consider using Locator
965    instead of manually manipulating process and point IDs.
966      */
967    ESCRIPT_DLL_API    ESCRIPT_DLL_API
968    void    const boost::python::tuple
969    calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;    maxGlobalDataPoint() const;
970    
971    
972    
973    /**    /**
974       \brief       \brief
975       Return the sign of each data point of this Data object.       Return the sign of each data point of this Data object.
# Line 882  class Data { Line 1029  class Data {
1029    /**    /**
1030       \brief       \brief
1031       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.
1032       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
1033       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
1034       first non-zero entry is positive.       first non-zero entry is positive.
1035       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
1036       *       *
# Line 1087  class Data { Line 1234  class Data {
1234    /**    /**
1235       \brief       \brief
1236       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.
1237        
1238       \param right Input - the power to raise the object to.       \param right Input - the power to raise the object to.
1239       *       *
1240     */     */
# Line 1098  class Data { Line 1245  class Data {
1245    /**    /**
1246       \brief       \brief
1247       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.
1248        
1249       \param left Input - the bases       \param left Input - the bases
1250       *       *
1251     */     */
# Line 1134  class Data { Line 1281  class Data {
1281    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1282    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1283    
1284      ESCRIPT_DLL_API
1285      Data& operator=(const Data& other);
1286    
1287    /**    /**
1288       \brief       \brief
1289       Overloaded operator -=       Overloaded operator -=
# Line 1226  class Data { Line 1376  class Data {
1376    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1377    inline    inline
1378    void    void
1379    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1380    
1381    /**    /**
1382       \brief       \brief
# Line 1237  class Data { Line 1387  class Data {
1387    */    */
1388    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1389    Data    Data
1390    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1391    
1392    /**    /**
1393       \brief       \brief
# Line 1250  class Data { Line 1400  class Data {
1400    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1401    void    void
1402    setSlice(const Data& value,    setSlice(const Data& value,
1403             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1404    
1405    /**    /**
1406       \brief       \brief
1407       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.  
1408    */    */
1409    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1410    void    void
1411    archiveData(const std::string fileName);          print(void);
1412    
1413    /**    /**
1414       \brief       \brief
1415       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1416       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1417       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.  
1418    */    */
1419    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1420    void          int
1421    extractData(const std::string fileName,          get_MPIRank(void) const;
               const FunctionSpace& fspace);  
   
1422    
1423    /**    /**
1424       \brief       \brief
1425       print the data values to stdout. Used for debugging       return the MPI rank number of the local data
1426                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1427                     is returned
1428    */    */
1429    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1430    void          int
1431      print(void);          get_MPISize(void) const;
1432    
1433    /**    /**
1434       \brief       \brief
1435       return the MPI rank number of the local data       return the MPI rank number of the local data
1436           MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()                   MPI_COMM_WORLD is assumed and returned.
          is returned  
1437    */    */
1438    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1439      int          MPI_Comm
1440      get_MPIRank(void) const;          get_MPIComm(void) const;
1441    
1442    /**    /**
1443       \brief       \brief
1444       return the MPI rank number of the local data       return the object produced by the factory, which is a DataConstant or DataExpanded
1445           MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()      TODO Ownership of this object should be explained in doco.
          is returned  
1446    */    */
1447    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1448      int          DataAbstract*
1449      get_MPISize(void) const;          borrowData(void) const;
1450    
   /**  
      \brief  
      return the MPI rank number of the local data  
          MPI_COMM_WORLD is assumed and returned.  
   */  
1451    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1452      MPI_Comm          DataAbstract_ptr
1453      get_MPIComm(void) const;          borrowDataPtr(void) const;
1454    
1455      ESCRIPT_DLL_API
1456            DataReady_ptr
1457            borrowReadyPtr(void) const;
1458    
1459    
1460    
1461    /**    /**
1462       \brief       \brief
1463       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.
1464         TODO Eventually these should be inlined.
1465         \param i - position(offset) in the underlying datastructure
1466    */    */
1467    
1468      ESCRIPT_DLL_API
1469            DataTypes::ValueType::const_reference
1470            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1471    
1472    
1473    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1474      DataAbstract*          DataTypes::ValueType::reference
1475      borrowData(void) const;          getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1476    
1477    
1478    
1479    /**
1480       \brief Create a buffer for use by getSample
1481       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1482       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1483      
1484       In multi-threaded sections, this needs to be called on each thread.
1485    
1486       \return A BufferGroup* if Data is lazy, NULL otherwise.
1487       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1488    */
1489      ESCRIPT_DLL_API
1490      BufferGroup*
1491      allocSampleBuffer() const;
1492    
1493    /**
1494       \brief Free a buffer allocated with allocSampleBuffer.
1495       \param buffer Input - pointer to the buffer to deallocate.
1496    */
1497    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1498    
1499   protected:   protected:
1500    
1501   private:   private:
1502    
1503      double
1504      LsupWorker() const;
1505    
1506      double
1507      supWorker() const;
1508    
1509      double
1510      infWorker() const;
1511    
1512      boost::python::object
1513      integrateWorker() const;
1514    
1515      void
1516      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1517    
1518      void
1519      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1520    
1521    
1522    /**    /**
1523       \brief       \brief
1524       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1395  class Data { Line 1590  class Data {
1590       \brief       \brief
1591       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1592    */    */
1593    template <class IValueType>  
1594    void    void
1595    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1596             const DataTypes::ShapeType& shape,
1597               const FunctionSpace& what,               const FunctionSpace& what,
1598               bool expanded);               bool expanded);
1599    
1600      void
1601      initialise(const WrappedArray& value,
1602                     const FunctionSpace& what,
1603                     bool expanded);
1604    
1605    //    //
1606    // flag to protect the data object against any update    // flag to protect the data object against any update
1607    bool m_protected;    bool m_protected;
1608      mutable bool m_shared;
1609      bool m_lazy;
1610    
1611    //    //
1612    // pointer to the actual data object    // pointer to the actual data object
1613    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1614      DataAbstract_ptr m_data;
1615    
1616    // If possible please use getReadyPtr instead.
1617    // But see warning below.
1618      const DataReady*
1619      getReady() const;
1620    
1621      DataReady*
1622      getReady();
1623    
1624    
1625    // Be wary of using this for local operations since it (temporarily) increases reference count.
1626    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1627    // getReady() instead
1628      DataReady_ptr
1629      getReadyPtr();
1630    
1631      const_DataReady_ptr
1632      getReadyPtr() const;
1633    
1634    
1635      /**
1636       \brief Update the Data's shared flag
1637       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1638       For internal use only.
1639      */
1640      void updateShareStatus(bool nowshared) const
1641      {
1642        m_shared=nowshared;     // m_shared is mutable
1643      }
1644    
1645      // In the isShared() method below:
1646      // A problem would occur if m_data (the address pointed to) were being modified
1647      // while the call m_data->is_shared is being executed.
1648      //
1649      // Q: So why do I think this code can be thread safe/correct?
1650      // A: We need to make some assumptions.
1651      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1652      // 2. We assume that no constructions or assignments which will share previously unshared
1653      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1654    //    //
1655    // pointer to the internal profiling data    // This means that the only transition we need to consider, is when a previously shared object is
1656    struct profDataEntry *profData;    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1657      // In those cases the m_shared flag changes to false after m_data has completed changing.
1658      // For any threads executing before the flag switches they will assume the object is still shared.
1659      bool isShared() const
1660      {
1661        return m_shared;
1662    /*  if (m_shared) return true;
1663        if (m_data->isShared())        
1664        {                  
1665            updateShareStatus(true);
1666            return true;
1667        }
1668        return false;*/
1669      }
1670    
1671      void forceResolve()
1672      {
1673        if (isLazy())
1674        {
1675            #ifdef _OPENMP
1676            if (omp_in_parallel())
1677            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1678            throw DataException("Please do not call forceResolve() in a parallel region.");
1679            }
1680            #endif
1681            resolve();
1682        }
1683      }
1684    
1685      /**
1686      \brief if another object is sharing out member data make a copy to work with instead.
1687      This code should only be called from single threaded sections of code.
1688      */
1689      void exclusiveWrite()
1690      {
1691    #ifdef _OPENMP
1692      if (omp_in_parallel())
1693      {
1694    // *((int*)0)=17;
1695        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1696      }
1697    #endif
1698        forceResolve();
1699        if (isShared())
1700        {
1701            DataAbstract* t=m_data->deepCopy();
1702            set_m_data(DataAbstract_ptr(t));
1703        }
1704      }
1705    
1706      /**
1707      \brief checks if caller can have exclusive write to the object
1708      */
1709      void checkExclusiveWrite()
1710      {
1711        if  (isLazy() || isShared())
1712        {
1713            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1714        }
1715      }
1716    
1717      /**
1718      \brief Modify the data abstract hosted by this Data object
1719      For internal use only.
1720      Passing a pointer to null is permitted (do this in the destructor)
1721      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1722      */
1723      void set_m_data(DataAbstract_ptr p);
1724    
1725      friend class DataAbstract;        // To allow calls to updateShareStatus
1726    
1727  };  };
1728    
1729  template <class IValueType>  }   // end namespace escript
1730  void  
1731  Data::initialise(const IValueType& value,  
1732                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1733                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1734    // so that I can dynamic cast between them below.
1735    #include "DataReady.h"
1736    #include "DataLazy.h"
1737    
1738    namespace escript
1739  {  {
1740    //  
1741    // Construct a Data object of the appropriate type.  inline
1742    // Construct the object first as there seems to be a bug which causes  const DataReady*
1743    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1744    // within the shared_ptr constructor.  {
1745    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1746      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1747      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
1748      m_data=temp_data;  }
1749    } else {  
1750      DataAbstract* temp=new DataConstant(value,what);  inline
1751      boost::shared_ptr<DataAbstract> temp_data(temp);  DataReady*
1752      m_data=temp_data;  Data::getReady()
1753    }  {
1754       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1755       EsysAssert((dr!=0), "Error - casting to DataReady.");
1756       return dr;
1757    }
1758    
1759    // Be wary of using this for local operations since it (temporarily) increases reference count.
1760    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1761    // getReady() instead
1762    inline
1763    DataReady_ptr
1764    Data::getReadyPtr()
1765    {
1766       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1767       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1768       return dr;
1769  }  }
1770    
1771    
1772    inline
1773    const_DataReady_ptr
1774    Data::getReadyPtr() const
1775    {
1776       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1777       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1778       return dr;
1779    }
1780    
1781    inline
1782    DataAbstract::ValueType::value_type*
1783    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1784    {
1785       if (isLazy())
1786       {
1787        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1788       }
1789       return getReady()->getSampleDataRW(sampleNo);
1790    }
1791    
1792    inline
1793    const DataAbstract::ValueType::value_type*
1794    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1795    {
1796       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1797       if (l!=0)
1798       {
1799        size_t offset=0;
1800        if (bufferg==NULL)
1801        {
1802            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1803        }
1804        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1805        return &((*res)[offset]);
1806       }
1807       return getReady()->getSampleDataRO(sampleNo);
1808    }
1809    
1810    
1811    
1812    /**
1813       Modify a filename for MPI parallel output to multiple files
1814    */
1815    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1816    
1817  /**  /**
1818     Binary Data object operators.     Binary Data object operators.
1819  */  */
1820  inline double rpow(double x,double y)  inline double rpow(double x,double y)
1821  {  {
1822      return pow(y,x);      return pow(y,x);
1823  };  }
1824    
1825  /**  /**
1826    \brief    \brief
# Line 1537  ESCRIPT_DLL_API Data operator*(const boo Line 1914  ESCRIPT_DLL_API Data operator*(const boo
1914  */  */
1915  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);
1916    
1917    
1918    
1919  /**  /**
1920    \brief    \brief
1921    Output operator    Output operator
# Line 1546  ESCRIPT_DLL_API std::ostream& operator<< Line 1925  ESCRIPT_DLL_API std::ostream& operator<<
1925  /**  /**
1926    \brief    \brief
1927    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1928    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1929    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1930    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1931    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1932  */  */
1933  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1934  Data  Data
1935  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1936                       Data& arg1,                       Data& arg_1,
1937                       int axis_offset=0,                       int axis_offset=0,
1938                       int transpose=0);                       int transpose=0);
1939    
1940  /**  /**
1941    \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  
1942    Perform the given binary operation with this and right as operands.    Perform the given binary operation with this and right as operands.
1943    Right is a Data object.    Right is a Data object.
1944  */  */
# Line 1579  Data::binaryOp(const Data& right, Line 1950  Data::binaryOp(const Data& right,
1950  {  {
1951     //     //
1952     // 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
1953     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1954       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.");
1955     }     }
1956    
1957       if (isLazy() || right.isLazy())
1958       {
1959         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1960       }
1961     //     //
1962     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1963     Data tempRight(right);     Data tempRight(right);
1964    
1965     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1966       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1967         //         //
1968         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1969         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1970       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1971         //         //
1972         // interpolate onto the RHS function space         // interpolate onto the RHS function space
1973         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
1974         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
1975           set_m_data(tempLeft.m_data);
1976       }       }
1977     }     }
1978     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1610  Data::binaryOp(const Data& right, Line 1988  Data::binaryOp(const Data& right,
1988       // of any data type       // of any data type
1989       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
1990       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
1991       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
1992     } else if (isTagged()) {     } else if (isTagged()) {
1993       //       //
1994       // 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 1632  Data::binaryOp(const Data& right, Line 2010  Data::binaryOp(const Data& right,
2010       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2011       escript::binaryOp(*leftC,*rightC,operation);       escript::binaryOp(*leftC,*rightC,operation);
2012     }     }
    #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);  
   }  
2013  }  }
2014    
2015  /**  /**
# Line 1707  Data::algorithm(BinaryFunction operation Line 2036  Data::algorithm(BinaryFunction operation
2036      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2037      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2038      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2039      } else if (isEmpty()) {
2040        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2041      } else if (isLazy()) {
2042        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2043      } else {
2044        throw DataException("Error - Data encapsulates an unknown type.");
2045    }    }
   return 0;  
2046  }  }
2047    
2048  /**  /**
2049    \brief    \brief
2050    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.
2051    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2052    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
2053    rank 0 Data object.    rank 0 Data object.
2054    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2055  */  */
# Line 1724  inline Line 2058  inline
2058  Data  Data
2059  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2060  {  {
2061    if (isExpanded()) {    if (isEmpty()) {
2062      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2063      }
2064      else if (isExpanded()) {
2065        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2066      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2067      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2068      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2069      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2070      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2071      return result;      return result;
2072    } else if (isTagged()) {    }
2073      else if (isTagged()) {
2074      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());  
2075      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2076      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2077        defval[0]=0;
2078        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2079      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2080      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2081    } else if (isConstant()) {    }
2082      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2083        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2084      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2085      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2086      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2087      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2088      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2089      return result;      return result;
2090      } else if (isLazy()) {
2091        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2092      } else {
2093        throw DataException("Error - Data encapsulates an unknown type.");
2094      }
2095    }
2096    
2097    /**
2098      \brief
2099      Compute a tensor operation with two Data objects
2100      \param arg_0 - Input - Data object
2101      \param arg_1 - Input - Data object
2102      \param operation - Input - Binary op functor
2103    */
2104    template <typename BinaryFunction>
2105    inline
2106    Data
2107    C_TensorBinaryOperation(Data const &arg_0,
2108                            Data const &arg_1,
2109                            BinaryFunction operation)
2110    {
2111      if (arg_0.isEmpty() || arg_1.isEmpty())
2112      {
2113         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2114      }
2115      if (arg_0.isLazy() || arg_1.isLazy())
2116      {
2117         throw DataException("Error - Operations not permitted on lazy data.");
2118      }
2119      // Interpolate if necessary and find an appropriate function space
2120      Data arg_0_Z, arg_1_Z;
2121      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2122        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2123          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2124          arg_1_Z = Data(arg_1);
2125        }
2126        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2127          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2128          arg_0_Z =Data(arg_0);
2129        }
2130        else {
2131          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2132        }
2133      } else {
2134          arg_0_Z = Data(arg_0);
2135          arg_1_Z = Data(arg_1);
2136      }
2137      // Get rank and shape of inputs
2138      int rank0 = arg_0_Z.getDataPointRank();
2139      int rank1 = arg_1_Z.getDataPointRank();
2140      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2141      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2142      int size0 = arg_0_Z.getDataPointSize();
2143      int size1 = arg_1_Z.getDataPointSize();
2144      // Declare output Data object
2145      Data res;
2146    
2147      if (shape0 == shape1) {
2148        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2149          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2150          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2151          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2152          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2153    
2154          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2155        }
2156        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2157    
2158          // Prepare the DataConstant input
2159          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2160    
2161          // Borrow DataTagged input from Data object
2162          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2163    
2164          // Prepare a DataTagged output 2
2165          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2166          res.tag();
2167          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2168    
2169          // Prepare offset into DataConstant
2170          int offset_0 = tmp_0->getPointOffset(0,0);
2171          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2172    
2173          // Get the pointers to the actual data
2174          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2175          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2176    
2177          // Compute a result for the default
2178          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2179          // Compute a result for each tag
2180          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2181          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2182          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2183            tmp_2->addTag(i->first);
2184            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2185            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2186    
2187            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2188          }
2189    
2190        }
2191        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2192          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2193          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2194          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2195          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2196    
2197          int sampleNo_1,dataPointNo_1;
2198          int numSamples_1 = arg_1_Z.getNumSamples();
2199          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2200          int offset_0 = tmp_0->getPointOffset(0,0);
2201          res.requireWrite();
2202          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2203          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2204            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2205              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2206              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2207              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2208              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2209              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2210              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2211            }
2212          }
2213    
2214        }
2215        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2216          // Borrow DataTagged input from Data object
2217          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2218    
2219          // Prepare the DataConstant input
2220          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2221    
2222          // Prepare a DataTagged output 2
2223          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2224          res.tag();
2225          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2226    
2227          // Prepare offset into DataConstant
2228          int offset_1 = tmp_1->getPointOffset(0,0);
2229    
2230          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2231          // Get the pointers to the actual data
2232          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2233          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2234          // Compute a result for the default
2235          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2236          // Compute a result for each tag
2237          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2238          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2239          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2240            tmp_2->addTag(i->first);
2241            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2242            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2243            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2244          }
2245    
2246        }
2247        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2248          // Borrow DataTagged input from Data object
2249          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2250    
2251          // Borrow DataTagged input from Data object
2252          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2253    
2254          // Prepare a DataTagged output 2
2255          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2256          res.tag();        // DataTagged output
2257          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2258    
2259          // Get the pointers to the actual data
2260          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2261          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2262          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2263    
2264          // Compute a result for the default
2265          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2266          // Merge the tags
2267          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2268          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2269          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2270          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2271            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2272          }
2273          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2274            tmp_2->addTag(i->first);
2275          }
2276          // Compute a result for each tag
2277          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2278          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2279    
2280            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2281            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2282            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2283    
2284            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2285          }
2286    
2287        }
2288        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2289          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2290          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2291          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2292          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2293          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2294    
2295          int sampleNo_0,dataPointNo_0;
2296          int numSamples_0 = arg_0_Z.getNumSamples();
2297          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2298          res.requireWrite();
2299          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2300          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2301            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2302            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2303            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2304              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2305              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2306              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2307              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2308              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2309            }
2310          }
2311    
2312        }
2313        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2314          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2315          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2316          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2317          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2318    
2319          int sampleNo_0,dataPointNo_0;
2320          int numSamples_0 = arg_0_Z.getNumSamples();
2321          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2322          int offset_1 = tmp_1->getPointOffset(0,0);
2323          res.requireWrite();
2324          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2325          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2326            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2327              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2328              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2329    
2330              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2331              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2332              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2333    
2334    
2335              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2336            }
2337          }
2338    
2339        }
2340        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2341          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2342          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2343          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2344          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2345          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2346    
2347          int sampleNo_0,dataPointNo_0;
2348          int numSamples_0 = arg_0_Z.getNumSamples();
2349          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2350          res.requireWrite();
2351          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2352          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2353            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2354            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2355            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2356              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2357              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2358              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2359              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2360              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2361            }
2362          }
2363    
2364        }
2365        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2366          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2367          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2368          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2369          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2370          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2371    
2372          int sampleNo_0,dataPointNo_0;
2373          int numSamples_0 = arg_0_Z.getNumSamples();
2374          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2375          res.requireWrite();
2376          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2377          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2378            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2379              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2380              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2381              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2382              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2383              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2384              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2385              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2386            }
2387          }
2388    
2389        }
2390        else {
2391          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2392        }
2393    
2394      } else if (0 == rank0) {
2395        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2396          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2397          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2398          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2399          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2400          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2401        }
2402        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2403    
2404          // Prepare the DataConstant input
2405          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2406    
2407          // Borrow DataTagged input from Data object
2408          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2409    
2410          // Prepare a DataTagged output 2
2411          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2412          res.tag();
2413          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2414    
2415          // Prepare offset into DataConstant
2416          int offset_0 = tmp_0->getPointOffset(0,0);
2417          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2418    
2419          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2420          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2421    
2422          // Compute a result for the default
2423          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2424          // Compute a result for each tag
2425          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2426          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2427          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2428            tmp_2->addTag(i->first);
2429            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2430            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2431            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2432          }
2433    
2434        }
2435        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2436    
2437          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2438          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2439          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2440          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2441    
2442          int sampleNo_1,dataPointNo_1;
2443          int numSamples_1 = arg_1_Z.getNumSamples();
2444          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2445          int offset_0 = tmp_0->getPointOffset(0,0);
2446          res.requireWrite();
2447          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2448          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2449            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2450              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2451              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2452              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2453              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2454              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2455              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2456    
2457            }
2458          }
2459    
2460        }
2461        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2462    
2463          // Borrow DataTagged input from Data object
2464          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2465    
2466          // Prepare the DataConstant input
2467          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2468    
2469          // Prepare a DataTagged output 2
2470          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2471          res.tag();
2472          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2473    
2474          // Prepare offset into DataConstant
2475          int offset_1 = tmp_1->getPointOffset(0,0);
2476          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2477    
2478          // Get the pointers to the actual data
2479          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2480          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2481    
2482    
2483          // Compute a result for the default
2484          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2485          // Compute a result for each tag
2486          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2487          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2488          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2489            tmp_2->addTag(i->first);
2490            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2491            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2492    
2493            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2494          }
2495    
2496        }
2497        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2498    
2499          // Borrow DataTagged input from Data object
2500          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2501    
2502          // Borrow DataTagged input from Data object
2503          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2504    
2505          // Prepare a DataTagged output 2
2506          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2507          res.tag();        // DataTagged output
2508          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2509    
2510          // Get the pointers to the actual data
2511          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2512          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2513          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2514    
2515          // Compute a result for the default
2516          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2517          // Merge the tags
2518          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2519          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2520          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2521          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2522            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2523          }
2524          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2525            tmp_2->addTag(i->first);
2526          }
2527          // Compute a result for each tag
2528          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2529          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2530            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2531            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2532            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2533    
2534            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2535          }
2536    
2537        }
2538        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2539    
2540          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2541          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2542          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2543          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2544          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2545    
2546          int sampleNo_0,dataPointNo_0;
2547          int numSamples_0 = arg_0_Z.getNumSamples();
2548          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2549          res.requireWrite();
2550          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2551          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2552            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2553            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2554            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2555              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2556              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2557              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2558              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2559              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2560            }
2561          }
2562    
2563        }
2564        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2565          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2566          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2567          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2568          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2569    
2570          int sampleNo_0,dataPointNo_0;
2571          int numSamples_0 = arg_0_Z.getNumSamples();
2572          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2573          int offset_1 = tmp_1->getPointOffset(0,0);
2574          res.requireWrite();
2575          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2576          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2577            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2578              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2579              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2580              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2581              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2582              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2583              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2584            }
2585          }
2586    
2587    
2588        }
2589        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2590    
2591          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2592          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2593          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2594          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2595          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2596    
2597          int sampleNo_0,dataPointNo_0;
2598          int numSamples_0 = arg_0_Z.getNumSamples();
2599          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2600          res.requireWrite();
2601          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2602          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2603            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2604            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2605            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2606              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2607              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2608              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2609              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2610              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2611            }
2612          }
2613    
2614        }
2615        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2616    
2617          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2618          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2619          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2620          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2621          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2622    
2623          int sampleNo_0,dataPointNo_0;
2624          int numSamples_0 = arg_0_Z.getNumSamples();
2625          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2626          res.requireWrite();
2627          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2628          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2629            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2630              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2631              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2632              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2633              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2634              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2635              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2636              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2637            }
2638          }
2639    
2640        }
2641        else {
2642          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2643        }
2644    
2645      } else if (0 == rank1) {
2646        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2647          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2648          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2649          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2650          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2651          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2652        }
2653        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2654    
2655          // Prepare the DataConstant input
2656          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2657    
2658          // Borrow DataTagged input from Data object
2659          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2660    
2661          // Prepare a DataTagged output 2
2662          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2663          res.tag();
2664          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2665    
2666          // Prepare offset into DataConstant
2667          int offset_0 = tmp_0->getPointOffset(0,0);
2668          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2669    
2670          //Get the pointers to the actual data
2671          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2672          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2673    
2674          // Compute a result for the default
2675          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2676          // Compute a result for each tag
2677          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2678          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2679          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2680            tmp_2->addTag(i->first);
2681            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2682            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2683            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2684          }
2685        }
2686        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2687    
2688          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2689          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2690          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2691          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2692    
2693          int sampleNo_1,dataPointNo_1;
2694          int numSamples_1 = arg_1_Z.getNumSamples();
2695          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2696          int offset_0 = tmp_0->getPointOffset(0,0);
2697          res.requireWrite();
2698          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2699          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2700            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2701              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2702              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2703              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2704              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2705              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2706              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2707            }
2708          }
2709    
2710        }
2711        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2712    
2713          // Borrow DataTagged input from Data object
2714          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2715    
2716          // Prepare the DataConstant input
2717          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2718    
2719          // Prepare a DataTagged output 2
2720          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2721          res.tag();
2722          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2723    
2724          // Prepare offset into DataConstant
2725          int offset_1 = tmp_1->getPointOffset(0,0);
2726          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2727          // Get the pointers to the actual data
2728          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2729          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2730          // Compute a result for the default
2731          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2732          // Compute a result for each tag
2733          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2734          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2735          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2736            tmp_2->addTag(i->first);
2737            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2738            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2739            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2740          }
2741    
2742        }
2743        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2744    
2745          // Borrow DataTagged input from Data object
2746          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2747    
2748          // Borrow DataTagged input from Data object
2749          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2750    
2751          // Prepare a DataTagged output 2
2752          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2753          res.tag();        // DataTagged output
2754          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2755    
2756          // Get the pointers to the actual data
2757          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2758          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2759          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2760    
2761          // Compute a result for the default
2762          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2763          // Merge the tags
2764          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2765          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2766          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2767          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2768            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2769          }
2770          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2771            tmp_2->addTag(i->first);
2772          }
2773          // Compute a result for each tag
2774          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2775          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2776            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2777            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2778            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2779            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2780          }
2781    
2782        }
2783        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2784    
2785          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2786          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2787          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2788          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2789          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2790    
2791          int sampleNo_0,dataPointNo_0;
2792          int numSamples_0 = arg_0_Z.getNumSamples();
2793          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2794          res.requireWrite();
2795          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2796          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2797            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2798            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2799            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2800              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2801              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2802              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2803              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2804              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2805            }
2806          }
2807    
2808        }
2809        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2810          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2811          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2812          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2813          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2814    
2815          int sampleNo_0,dataPointNo_0;
2816          int numSamples_0 = arg_0_Z.getNumSamples();
2817          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2818          int offset_1 = tmp_1->getPointOffset(0,0);
2819          res.requireWrite();
2820          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2821          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2822            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2823              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2824              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2825              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2826              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2827              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2828              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2829            }
2830          }
2831    
2832    
2833        }
2834        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2835    
2836          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2837          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2838          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2839          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2840          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2841    
2842          int sampleNo_0,dataPointNo_0;
2843          int numSamples_0 = arg_0_Z.getNumSamples();
2844          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2845          res.requireWrite();
2846          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2847          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2848            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2849            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2850            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2851              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2852              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2853              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2854              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2855              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2856            }
2857          }
2858    
2859        }
2860        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2861    
2862          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2863          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2864          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2865          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2866          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2867    
2868          int sampleNo_0,dataPointNo_0;
2869          int numSamples_0 = arg_0_Z.getNumSamples();
2870          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2871          res.requireWrite();
2872          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2873          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2874            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2875              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2876              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2877              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2878              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2879              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2880              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2881              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2882            }
2883          }
2884    
2885        }
2886        else {
2887          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2888        }
2889    
2890      } else {
2891        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2892      }
2893    
2894      return res;
2895    }
2896    
2897    template <typename UnaryFunction>
2898    Data
2899    C_TensorUnaryOperation(Data const &arg_0,
2900                           UnaryFunction operation)
2901    {
2902      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2903      {
2904         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2905      }
2906      if (arg_0.isLazy())
2907      {
2908         throw DataException("Error - Operations not permitted on lazy data.");
2909    }    }
2910    Data falseRetVal; // to keep compiler quiet    // Interpolate if necessary and find an appropriate function space
2911    return falseRetVal;    Data arg_0_Z = Data(arg_0);
2912    
2913      // Get rank and shape of inputs
2914      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2915      int size0 = arg_0_Z.getDataPointSize();
2916    
2917      // Declare output Data object
2918      Data res;
2919    
2920      if (arg_0_Z.isConstant()) {
2921        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2922        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2923        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2924        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2925      }
2926      else if (arg_0_Z.isTagged()) {
2927    
2928        // Borrow DataTagged input from Data object
2929        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2930    
2931        // Prepare a DataTagged output 2
2932        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2933        res.tag();
2934        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2935    
2936        // Get the pointers to the actual data
2937        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2938        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2939        // Compute a result for the default
2940        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2941        // Compute a result for each tag
2942        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2943        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2944        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2945          tmp_2->addTag(i->first);
2946          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2947          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2948          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2949        }
2950    
2951      }
2952      else if (arg_0_Z.isExpanded()) {
2953    
2954        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2955        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2956        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2957    
2958        int sampleNo_0,dataPointNo_0;
2959        int numSamples_0 = arg_0_Z.getNumSamples();
2960        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2961        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2962        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2963          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2964            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2965            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2966            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2967            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2968            tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2969          }
2970        }
2971      }
2972      else {
2973        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
2974      }
2975    
2976      return res;
2977  }  }
2978    
2979  }  }

Legend:
Removed from v.922  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26