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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26