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

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

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

revision 854 by gross, Thu Sep 21 05:29:42 2006 UTC revision 2635 by jfenwick, Thu Aug 27 04:54:41 2009 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2009 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
8   *  Licensed under the Open Software License version 3.0    *  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Licensed under the Open Software License version 3.0
10   *                                                          *  * http://www.opensource.org/licenses/osl-3.0.php
11   ************************************************************  *
12  */  *******************************************************/
13    
14    
15  /** \file Data.h */  /** \file Data.h */
16    
# Line 17  Line 18 
18  #define DATA_H  #define DATA_H
19  #include "system_dep.h"  #include "system_dep.h"
20    
21    #include "DataTypes.h"
22  #include "DataAbstract.h"  #include "DataAbstract.h"
23  #include "DataAlgorithm.h"  #include "DataAlgorithm.h"
24  #include "FunctionSpace.h"  #include "FunctionSpace.h"
# Line 24  Line 26 
26  #include "UnaryOp.h"  #include "UnaryOp.h"
27  #include "DataException.h"  #include "DataException.h"
28    
29    
30  extern "C" {  extern "C" {
31  #include "DataC.h"  #include "DataC.h"
32  #include "paso/Paso.h"  //#include <omp.h>
33  }  }
34    
35  #ifndef PASO_MPI  #include "esysmpi.h"
 #define MPI_Comm long  
 #endif  
   
36  #include <string>  #include <string>
37  #include <algorithm>  #include <algorithm>
38    #include <sstream>
39    
40  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
41  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
42  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
43  #include <boost/python/numeric.hpp>  
44    #include "BufferGroup.h"
45    
46  namespace escript {  namespace escript {
47    
# Line 48  namespace escript { Line 50  namespace escript {
50  class DataConstant;  class DataConstant;
51  class DataTagged;  class DataTagged;
52  class DataExpanded;  class DataExpanded;
53    class DataLazy;
54    
55  /**  /**
56     \brief     \brief
57     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
58    
59     Description:     Description:
60     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
61     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
62     the object created for the lifetime of the object.     of the Data object.
63     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
64     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
65       Doing so will lead to invalid memory access.
66       This should not affect any methods exposed via boost::python.
67  */  */
68  class Data {  class Data {
69    
# Line 102  class Data { Line 106  class Data {
106         const FunctionSpace& what);         const FunctionSpace& what);
107    
108    /**    /**
109       \brief      \brief Copy Data from an existing vector
110       Constructor which copies data from a DataArrayView.    */
111    
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
112    ESCRIPT_DLL_API    ESCRIPT_DLL_API
113    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
114         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
115         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
116                     bool expanded=false);
117    
118    /**    /**
119       \brief       \brief
120       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
121    
122       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
123       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 129  class Data { Line 128  class Data {
128    */    */
129    ESCRIPT_DLL_API    ESCRIPT_DLL_API
130    Data(double value,    Data(double value,
131         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
132         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
133         bool expanded=false);         bool expanded=false);
134    
# Line 142  class Data { Line 141  class Data {
141    */    */
142    ESCRIPT_DLL_API    ESCRIPT_DLL_API
143    Data(const Data& inData,    Data(const Data& inData,
144         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
145    
146    /**    /**
147       \brief       \brief
148       Constructor which will create Tagged data if expanded is false.       Constructor which copies data from any object that can be treated like a python array/sequence.
      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);  
   
   /**  
      \brief  
      Constructor which copies data from any object that can be converted into  
      a python numarray.  
149    
150       \param value - Input - Input data.       \param value - Input - Input data.
151       \param what - Input - A description of what this data represents.       \param what - Input - A description of what this data represents.
# Line 199  class Data { Line 161  class Data {
161    /**    /**
162       \brief       \brief
163       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
164       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
165       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
166    
167       \param value - Input - Input data.       \param value - Input - Input data.
168       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
# Line 214  class Data { Line 176  class Data {
176       Constructor which creates a DataConstant of "shape" with constant value.       Constructor which creates a DataConstant of "shape" with constant value.
177    */    */
178    ESCRIPT_DLL_API    ESCRIPT_DLL_API
179    Data(double value,    Data(double value,
180         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
181         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
182         bool expanded=false);         bool expanded=false);
183    
184    
185    
186      /**
187        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
188        Once you have passed the pointer, do not delete it.
189      */
190      ESCRIPT_DLL_API
191      explicit Data(DataAbstract* underlyingdata);
192    
193      /**
194        \brief Create a Data based on the supplied DataAbstract
195      */
196      ESCRIPT_DLL_API
197      explicit Data(DataAbstract_ptr underlyingdata);
198    
199    /**    /**
200       \brief       \brief
201       Destructor       Destructor
# Line 226  class Data { Line 204  class Data {
204    ~Data();    ~Data();
205    
206    /**    /**
207       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
208    */    */
209    ESCRIPT_DLL_API    ESCRIPT_DLL_API
210    void    void
211    copy(const Data& other);    copy(const Data& other);
212    
213    /**    /**
214         \brief Return a pointer to a deep copy of this object.
215      */
216      ESCRIPT_DLL_API
217      Data
218      copySelf();
219    
220    
221      /**
222         \brief produce a delayed evaluation version of this Data.
223      */
224      ESCRIPT_DLL_API
225      Data
226      delay();
227    
228      /**
229         \brief convert the current data into lazy data.
230      */
231      ESCRIPT_DLL_API
232      void
233      delaySelf();
234    
235    
236      /**
237       Member access methods.       Member access methods.
238    */    */
239    
240    /**    /**
241       \brief       \brief
242       switches on update protection       switches on update protection
243    
244    */    */
245    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 248  class Data { Line 248  class Data {
248    
249    /**    /**
250       \brief       \brief
251       Returns trueif the data object is protected against update       Returns true, if the data object is protected against update
252    
253    */    */
254    ESCRIPT_DLL_API    ESCRIPT_DLL_API
255    bool    bool
256    isProtected() const;    isProtected() const;
257    
258    
259    /**
260       \brief
261       Return the value of a data point as a python tuple.
262    */
263      ESCRIPT_DLL_API
264      const boost::python::object
265      getValueOfDataPointAsTuple(int dataPointNo);
266    
267    /**    /**
268       \brief       \brief
269       Return the values of all data-points as a single python numarray object.       sets the values of a data-point from a python object on this process
270    */    */
271    ESCRIPT_DLL_API    ESCRIPT_DLL_API
272    const boost::python::numeric::array    void
273    convertToNumArray();    setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
274    
275    /**    /**
276       \brief       \brief
277       Return the values of all data-points for the given sample as a single python numarray object.       sets the values of a data-point from a array-like object on this process
278    */    */
279    ESCRIPT_DLL_API    ESCRIPT_DLL_API
280    const boost::python::numeric::array    void
281    convertToNumArrayFromSampleNo(int sampleNo);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282    
283    /**    /**
284       \brief       \brief
285       Return the value of the specified data-point as a single python numarray object.       sets the values of a data-point on this process
286    */    */
 #ifndef PASO_MPI    
287    ESCRIPT_DLL_API    ESCRIPT_DLL_API
288    const boost::python::numeric::array    void
289    convertToNumArrayFromDPNo(int ProcNo,    setValueOfDataPoint(int dataPointNo, const double);
                                                         int sampleNo,  
                             int dataPointNo);  
 #else  
   ESCRIPT_DLL_API  
   const boost::python::numeric::array  
   convertToNumArrayFromDPNo(int procNo,  
                 int sampleNo,  
                 int dataPointNo);  
 #endif  
   
290    
291    /**    /**
292       \brief       \brief Return a data point across all processors as a python tuple.
      Fills the expanded Data object from values of a python numarray object.  
293    */    */
294    ESCRIPT_DLL_API    ESCRIPT_DLL_API
295    void    const boost::python::object
296    fillFromNumArray(const boost::python::numeric::array);    getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
297    
298    /**    /**
299       \brief       \brief
300       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
301    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
302    */    */
303    ESCRIPT_DLL_API    ESCRIPT_DLL_API
304    int    int
# Line 316  class Data { Line 312  class Data {
312    escriptDataC    escriptDataC
313    getDataC();    getDataC();
314    
315    
316    
317    /**    /**
318       \brief       \brief
319       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 325  class Data { Line 323  class Data {
323    getDataC() const;    getDataC() const;
324    
325    /**    /**
326       \brief      \brief How much space is required to evaulate a sample of the Data.
      Write the data as a string.  
327    */    */
328    ESCRIPT_DLL_API    ESCRIPT_DLL_API
329    inline    size_t
330    std::string    getSampleBufferSize() const;
331    toString() const  
332    {  
     return m_data->toString();  
   }  
333    
334    /**    /**
335       \brief       \brief
336       Return the DataArrayView of the point data. This essentially contains       Write the data as a string. For large amounts of data, a summary is printed.
      the shape information for each data point although it also may be used  
      to manipulate the point data.  
337    */    */
338    ESCRIPT_DLL_API    ESCRIPT_DLL_API
339    inline    std::string
340    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
341    
342    /**    /**
343       \brief       \brief
# Line 363  class Data { Line 352  class Data {
352       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
353       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
354       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
355    */    */
356    ESCRIPT_DLL_API    ESCRIPT_DLL_API
357    void    void
358    tag();    tag();
359    
360    /**    /**
361        \brief If this data is lazy, then convert it to ready data.
362        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
363      */
364      ESCRIPT_DLL_API
365      void
366      resolve();
367    
368    
369      /**
370       \brief Ensures data is ready for write access.
371      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
372      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
373      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
374      Doing so might introduce additional sharing.
375      */
376      ESCRIPT_DLL_API
377      void
378      requireWrite();
379    
380      /**
381       \brief       \brief
382       Return true if this Data is expanded.       Return true if this Data is expanded.
383         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
384    */    */
385    ESCRIPT_DLL_API    ESCRIPT_DLL_API
386    bool    bool
# Line 379  class Data { Line 388  class Data {
388    
389    /**    /**
390       \brief       \brief
391         Return true if this Data is expanded or resolves to expanded.
392         That is, if it has a separate value for each datapoint in the sample.
393      */
394      ESCRIPT_DLL_API
395      bool
396      actsExpanded() const;
397      
398    
399      /**
400         \brief
401       Return true if this Data is tagged.       Return true if this Data is tagged.
402    */    */
403    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 395  class Data { Line 414  class Data {
414    
415    /**    /**
416       \brief       \brief
417       Return true if this Data is empty.       Return true if this Data is constant or resolves to constant.
418         That is, if there is only one value stored, regardless of how many
419         datapoints there are per sample.
420      */
421      ESCRIPT_DLL_API
422      bool
423      actsConstant() const;
424    
425      /**
426         \brief Return true if this Data is lazy.
427      */
428      ESCRIPT_DLL_API
429      bool
430      isLazy() const;
431    
432      /**
433         \brief Return true if this data is ready.
434      */
435      ESCRIPT_DLL_API
436      bool
437      isReady() const;
438    
439      /**
440         \brief
441         Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
442    contains datapoints.
443    */    */
444    ESCRIPT_DLL_API    ESCRIPT_DLL_API
445    bool    bool
# Line 427  class Data { Line 471  class Data {
471    */    */
472    ESCRIPT_DLL_API    ESCRIPT_DLL_API
473    inline    inline
474    const AbstractDomain&  //   const AbstractDomain&
475      const_Domain_ptr
476    getDomain() const    getDomain() const
477    {    {
478       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
479    }    }
480    
481    
482      /**
483         \brief
484         Return the domain.
485         TODO: For internal use only.   This should be removed.
486      */
487      ESCRIPT_DLL_API
488      inline
489    //   const AbstractDomain&
490      Domain_ptr
491      getDomainPython() const
492      {
493         return getFunctionSpace().getDomainPython();
494      }
495    
496    /**    /**
497       \brief       \brief
498       Return a copy of the domain.       Return a copy of the domain.
# Line 447  class Data { Line 507  class Data {
507    */    */
508    ESCRIPT_DLL_API    ESCRIPT_DLL_API
509    inline    inline
510    int    unsigned int
511    getDataPointRank() const    getDataPointRank() const
512    {    {
513      return m_data->getPointDataView().getRank();      return m_data->getRank();
514    }    }
515    
516    /**    /**
517       \brief       \brief
518         Return the number of data points
519      */
520      ESCRIPT_DLL_API
521      inline
522      int
523      getNumDataPoints() const
524      {
525        return getNumSamples() * getNumDataPointsPerSample();
526      }
527      /**
528         \brief
529       Return the number of samples.       Return the number of samples.
530    */    */
531    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 477  class Data { Line 548  class Data {
548      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
549    }    }
550    
551    
552      /**
553        \brief
554        Return the number of values in the shape for this object.
555      */
556      ESCRIPT_DLL_API
557      int
558      getNoValues() const
559      {
560        return m_data->getNoValues();
561      }
562    
563    
564      /**
565         \brief
566         dumps the object into a netCDF file
567      */
568      ESCRIPT_DLL_API
569      void
570      dump(const std::string fileName) const;
571    
572     /**
573      \brief returns the values of the object as a list of tuples (one for each datapoint).
574    
575      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
576    If false, the result is a list of scalars [1, 2, ...]
577     */
578      ESCRIPT_DLL_API
579      const boost::python::object
580      toListOfTuples(bool scalarastuple=true);
581    
582    
583     /**
584        \brief
585        Return the sample data for the given sample no. This is not the
586        preferred interface but is provided for use by C code.
587        The bufferg parameter is only required for LazyData.
588        \param sampleNo - Input - the given sample no.
589        \param bufferg - A buffer to compute (and store) sample data in will be selected from this group.
590        \return pointer to the sample data.
591    */
592      ESCRIPT_DLL_API
593      inline
594      const DataAbstract::ValueType::value_type*
595      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg=0);
596    
597    
598    /**    /**
599       \brief       \brief
600       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
601       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
602       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
603         \return pointer to the sample data.
604    */    */
605    ESCRIPT_DLL_API    ESCRIPT_DLL_API
606    inline    inline
607    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
608    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
609    {  
     return m_data->getSampleData(sampleNo);  
   }  
610    
611    /**    /**
612       \brief       \brief
# Line 507  class Data { Line 624  class Data {
624    
625    /**    /**
626       \brief       \brief
627       Assign the given value to the data-points referenced by the given       Return a reference into the DataVector which points to the specified data point.
628       reference number.       \param sampleNo - Input -
629         \param dataPointNo - Input -
      The value supplied is a python numarray object.  The data from this numarray  
      is unpacked into a DataArray, and this is used to set the corresponding  
      data-points in the underlying Data object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Input - value to assign to data-points associated with  
                             the given reference number.  
630    */    */
631    ESCRIPT_DLL_API    ESCRIPT_DLL_API
632    void    DataTypes::ValueType::const_reference
633    setRefValue(int ref,    getDataPointRO(int sampleNo, int dataPointNo);
               const boost::python::numeric::array& value);  
634    
635    /**    /**
636       \brief       \brief
637       Return the values associated with the data-points referenced by the given       Return a reference into the DataVector which points to the specified data point.
638       reference number.       \param sampleNo - Input -
639         \param dataPointNo - Input -
      The value supplied is a python numarray object. The data from the corresponding  
      data-points in this Data object are packed into the given numarray object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Output - object to receive values from data-points  
                              associated with the given reference number.  
640    */    */
641    ESCRIPT_DLL_API    ESCRIPT_DLL_API
642    void    DataTypes::ValueType::reference
643    getRefValue(int ref,    getDataPointRW(int sampleNo, int dataPointNo);
644                boost::python::numeric::array& value);  
645    
646    
647    /**    /**
648       \brief       \brief
649       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 -  
650    */    */
651    ESCRIPT_DLL_API    ESCRIPT_DLL_API
652    inline    inline
653    DataArrayView    DataTypes::ValueType::size_type
654    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
655                 int dataPointNo)                 int dataPointNo)
656    {    {
657          return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
658    }    }
659    
660    /**    /**
# Line 568  class Data { Line 662  class Data {
662       Return a reference to the data point shape.       Return a reference to the data point shape.
663    */    */
664    ESCRIPT_DLL_API    ESCRIPT_DLL_API
665    const DataArrayView::ShapeType&    inline
666    getDataPointShape() const;    const DataTypes::ShapeType&
667      getDataPointShape() const
668      {
669        return m_data->getShape();
670      }
671    
672    /**    /**
673       \brief       \brief
# Line 593  class Data { Line 691  class Data {
691       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
692    */    */
693    ESCRIPT_DLL_API    ESCRIPT_DLL_API
694    DataArrayView::ValueType::size_type    DataTypes::ValueType::size_type
695    getLength() const;    getLength() const;
696    
697    
698    
699    /**    /**
700       \brief       \brief
701       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
702       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
703       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
704         \param name - Input - name of tag.
705         \param value - Input - Value to associate with given key.
706      */
707      ESCRIPT_DLL_API
708      void
709      setTaggedValueByName(std::string name,
710                           const boost::python::object& value);
711    
712      /**
713         \brief
714         Assign the given value to the tag. Implicitly converts this
715         object to type DataTagged if it is constant.
716    
717       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
718       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
719      ==>*      ==>*
# Line 613  class Data { Line 726  class Data {
726    /**    /**
727       \brief       \brief
728       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
729       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
730       cannot be converted to a DataTagged object.  
731       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
732         \param pointshape - Input - The shape of the value parameter
733       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
734      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
735    */    */
736    ESCRIPT_DLL_API    ESCRIPT_DLL_API
737    void    void
738    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
739                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
740                            const DataTypes::ValueType& value,
741                int dataOffset=0);
742    
743    
744    
745    /**    /**
746      \brief      \brief
# Line 639  class Data { Line 757  class Data {
757    
758    /**    /**
759       \brief       \brief
760         set all values to zero
761         *
762      */
763      ESCRIPT_DLL_API
764      void
765      setToZero();
766    
767      /**
768         \brief
769       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
770       the result as a Data object.       the result as a Data object.
771       *       *
# Line 647  class Data { Line 774  class Data {
774    Data    Data
775    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
776    
777    
778      ESCRIPT_DLL_API
779      Data
780      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
781                           double undef, Data& B, double Bmin, double Bstep);
782    
783    
784      // This signature needs to be tuned
785      ESCRIPT_DLL_API
786      Data
787      interpolateFromTable(boost::python::object table, double Amin, double Astep,
788                           double undef, Data& B, double Bmin, double Bstep);
789    
790    /**    /**
791       \brief       \brief
792       Calculates the gradient of the data at the data points of functionspace.       Calculates the gradient of the data at the data points of functionspace.
# Line 662  class Data { Line 802  class Data {
802    grad() const;    grad() const;
803    
804    /**    /**
805       \brief      \brief
806       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
807       *    */
808      ESCRIPT_DLL_API
809      boost::python::object
810      integrateToTuple_const() const;
811    
812    
813      /**
814        \brief
815         Calculate the integral over the function space domain as a python tuple.
816    */    */
817    ESCRIPT_DLL_API    ESCRIPT_DLL_API
818    boost::python::numeric::array    boost::python::object
819    integrate() const;    integrateToTuple();
820    
821    
822    
823    /**    /**
824       \brief       \brief
# Line 735  class Data { Line 885  class Data {
885    /**    /**
886       \brief       \brief
887       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
888       *  
889         The method is not const because lazy data needs to be expanded before Lsup can be computed.
890         The _const form can be used when the Data object is const, however this will only work for
891         Data which is not Lazy.
892    
893         For Data which contain no samples (or tagged Data for which no tags in use have a value)
894         zero is returned.
895    */    */
896    ESCRIPT_DLL_API    ESCRIPT_DLL_API
897    double    double
898    Lsup() const;    Lsup();
899    
   /**  
      \brief  
      Return the minimum absolute value of this Data object.  
      *  
   */  
900    ESCRIPT_DLL_API    ESCRIPT_DLL_API
901    double    double
902    Linf() const;    Lsup_const() const;
903    
904    
905    /**    /**
906       \brief       \brief
907       Return the maximum value of this Data object.       Return the maximum value of this Data object.
908       *  
909         The method is not const because lazy data needs to be expanded before sup can be computed.
910         The _const form can be used when the Data object is const, however this will only work for
911         Data which is not Lazy.
912    
913         For Data which contain no samples (or tagged Data for which no tags in use have a value)
914         a large negative value is returned.
915    */    */
916    ESCRIPT_DLL_API    ESCRIPT_DLL_API
917    double    double
918    sup() const;    sup();
919    
920      ESCRIPT_DLL_API
921      double
922      sup_const() const;
923    
924    
925    /**    /**
926       \brief       \brief
927       Return the minimum value of this Data object.       Return the minimum value of this Data object.
928       *  
929         The method is not const because lazy data needs to be expanded before inf can be computed.
930         The _const form can be used when the Data object is const, however this will only work for
931         Data which is not Lazy.
932    
933         For Data which contain no samples (or tagged Data for which no tags in use have a value)
934         a large positive value is returned.
935    */    */
936    ESCRIPT_DLL_API    ESCRIPT_DLL_API
937    double    double
938    inf() const;    inf();
939    
940      ESCRIPT_DLL_API
941      double
942      inf_const() const;
943    
944    
945    
946    /**    /**
947       \brief       \brief
# Line 798  class Data { Line 973  class Data {
973    /**    /**
974       \brief       \brief
975       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
976       the minimum value in this Data object.       the minimum component value in this Data object.
977         \note If you are working in python, please consider using Locator
978    instead of manually manipulating process and point IDs.
979    */    */
980    ESCRIPT_DLL_API    ESCRIPT_DLL_API
981    const boost::python::tuple    const boost::python::tuple
982    mindp() const;    minGlobalDataPoint() const;
983    
984      /**
985         \brief
986         Return the (sample number, data-point number) of the data point with
987         the minimum component value in this Data object.
988         \note If you are working in python, please consider using Locator
989    instead of manually manipulating process and point IDs.
990      */
991    ESCRIPT_DLL_API    ESCRIPT_DLL_API
992    void    const boost::python::tuple
993    calc_mindp(int& ProcNo,    maxGlobalDataPoint() const;
994                          int& SampleNo,    
995               int& DataPointNo) const;  
996    
997    /**    /**
998       \brief       \brief
999       Return the sign of each data point of this Data object.       Return the sign of each data point of this Data object.
# Line 868  class Data { Line 1053  class Data {
1053    /**    /**
1054       \brief       \brief
1055       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.
1056       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
1057       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
1058       first non-zero entry is positive.       first non-zero entry is positive.
1059       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
1060       *       *
# Line 889  class Data { Line 1074  class Data {
1074    
1075    /**    /**
1076       \brief       \brief
1077         Return the error function erf of each data point of this Data object.
1078         *
1079      */
1080      ESCRIPT_DLL_API
1081      Data
1082      erf() const;
1083    
1084      /**
1085         \brief
1086       Return the sin of each data point of this Data object.       Return the sin of each data point of this Data object.
1087       *       *
1088    */    */
# Line 1064  class Data { Line 1258  class Data {
1258    /**    /**
1259       \brief       \brief
1260       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.
1261        
1262       \param right Input - the power to raise the object to.       \param right Input - the power to raise the object to.
1263       *       *
1264     */     */
# Line 1075  class Data { Line 1269  class Data {
1269    /**    /**
1270       \brief       \brief
1271       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.
1272        
1273       \param left Input - the bases       \param left Input - the bases
1274       *       *
1275     */     */
# Line 1100  class Data { Line 1294  class Data {
1294    void    void
1295    saveVTK(std::string fileName) const;    saveVTK(std::string fileName) const;
1296    
1297    
1298    
1299    /**    /**
1300       \brief       \brief
1301       Overloaded operator +=       Overloaded operator +=
# Line 1111  class Data { Line 1307  class Data {
1307    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1308    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1309    
1310      ESCRIPT_DLL_API
1311      Data& operator=(const Data& other);
1312    
1313    /**    /**
1314       \brief       \brief
1315       Overloaded operator -=       Overloaded operator -=
# Line 1203  class Data { Line 1402  class Data {
1402    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1403    inline    inline
1404    void    void
1405    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1406    
1407    /**    /**
1408       \brief       \brief
# Line 1214  class Data { Line 1413  class Data {
1413    */    */
1414    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1415    Data    Data
1416    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1417    
1418    /**    /**
1419       \brief       \brief
# Line 1227  class Data { Line 1426  class Data {
1426    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1427    void    void
1428    setSlice(const Data& value,    setSlice(const Data& value,
1429             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1430    
1431    /**    /**
1432       \brief       \brief
1433       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.  
1434    */    */
1435    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1436    void    void
1437    archiveData(const std::string fileName);          print(void);
1438    
1439    /**    /**
1440       \brief       \brief
1441       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1442       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1443       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.  
1444    */    */
1445    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1446    void          int
1447    extractData(const std::string fileName,          get_MPIRank(void) const;
               const FunctionSpace& fspace);  
   
1448    
1449    /**    /**
1450       \brief       \brief
1451       print the data values to stdout. Used for debugging       return the MPI rank number of the local data
1452                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1453                     is returned
1454    */    */
1455    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1456    void          int
1457      print(void);          get_MPISize(void) const;
1458    
1459    /**    /**
1460       \brief       \brief
1461       return the MPI rank number of the local data       return the MPI rank number of the local data
1462           MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()                   MPI_COMM_WORLD is assumed and returned.
          is returned  
1463    */    */
1464    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1465      int          MPI_Comm
1466      get_MPIRank(void) const;          get_MPIComm(void) const;
1467    
1468    /**    /**
1469       \brief       \brief
1470       return the MPI rank number of the local data       return the object produced by the factory, which is a DataConstant or DataExpanded
1471           MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()      TODO Ownership of this object should be explained in doco.
          is returned  
1472    */    */
1473    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1474      int          DataAbstract*
1475      get_MPISize(void) const;          borrowData(void) const;
1476    
1477      ESCRIPT_DLL_API
1478            DataAbstract_ptr
1479            borrowDataPtr(void) const;
1480    
   /**  
      \brief  
      return the MPI rank number of the local data  
          MPI_COMM_WORLD is assumed and returned.  
   */  
1481    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1482      MPI_Comm          DataReady_ptr
1483      get_MPIComm(void) const;          borrowReadyPtr(void) const;
1484    
1485    
1486    
1487    /**    /**
1488       \brief       \brief
1489       return the object produced by the factory, which is a DataConstant or DataExpanded       Return a pointer to the beginning of the datapoint at the specified offset.
1490         TODO Eventually these should be inlined.
1491         \param i - position(offset) in the underlying datastructure
1492    */    */
1493    
1494      ESCRIPT_DLL_API
1495            DataTypes::ValueType::const_reference
1496            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1497    
1498    
1499    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1500      DataAbstract*          DataTypes::ValueType::reference
1501      borrowData(void) const;          getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1502    
1503    
1504    
1505    /**
1506       \brief Create a buffer for use by getSample
1507       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1508       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1509      
1510       In multi-threaded sections, this needs to be called on each thread.
1511    
1512       \return A BufferGroup* if Data is lazy, NULL otherwise.
1513       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1514    */
1515      ESCRIPT_DLL_API
1516      BufferGroup*
1517      allocSampleBuffer() const;
1518    
1519    /**
1520       \brief Free a buffer allocated with allocSampleBuffer.
1521       \param buffer Input - pointer to the buffer to deallocate.
1522    */
1523    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1524    
1525   protected:   protected:
1526    
1527   private:   private:
1528    
1529      double
1530      LsupWorker() const;
1531    
1532      double
1533      supWorker() const;
1534    
1535      double
1536      infWorker() const;
1537    
1538      boost::python::object
1539      integrateWorker() const;
1540    
1541      void
1542      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1543    
1544      void
1545      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1546    
1547    
1548    /**    /**
1549       \brief       \brief
1550       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1372  class Data { Line 1616  class Data {
1616       \brief       \brief
1617       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1618    */    */
1619    template <class IValueType>  
1620    void    void
1621    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1622             const DataTypes::ShapeType& shape,
1623               const FunctionSpace& what,               const FunctionSpace& what,
1624               bool expanded);               bool expanded);
1625    
1626      void
1627      initialise(const WrappedArray& value,
1628                     const FunctionSpace& what,
1629                     bool expanded);
1630    
1631    //    //
1632    // flag to protect the data object against any update    // flag to protect the data object against any update
1633    bool m_protected;    bool m_protected;
1634      mutable bool m_shared;
1635      bool m_lazy;
1636    
1637    //    //
1638    // pointer to the actual data object    // pointer to the actual data object
1639    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1640      DataAbstract_ptr m_data;
1641    
1642    // If possible please use getReadyPtr instead.
1643    // But see warning below.
1644      const DataReady*
1645      getReady() const;
1646    
1647      DataReady*
1648      getReady();
1649    
1650    
1651    // Be wary of using this for local operations since it (temporarily) increases reference count.
1652    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1653    // getReady() instead
1654      DataReady_ptr
1655      getReadyPtr();
1656    
1657      const_DataReady_ptr
1658      getReadyPtr() const;
1659    
1660    
1661      /**
1662       \brief Update the Data's shared flag
1663       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1664       For internal use only.
1665      */
1666      void updateShareStatus(bool nowshared) const
1667      {
1668        m_shared=nowshared;     // m_shared is mutable
1669      }
1670    
1671      // In the isShared() method below:
1672      // A problem would occur if m_data (the address pointed to) were being modified
1673      // while the call m_data->is_shared is being executed.
1674      //
1675      // Q: So why do I think this code can be thread safe/correct?
1676      // A: We need to make some assumptions.
1677      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1678      // 2. We assume that no constructions or assignments which will share previously unshared
1679      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1680    //    //
1681    // pointer to the internal profiling data    // This means that the only transition we need to consider, is when a previously shared object is
1682    struct profDataEntry *profData;    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1683      // In those cases the m_shared flag changes to false after m_data has completed changing.
1684      // For any threads executing before the flag switches they will assume the object is still shared.
1685      bool isShared() const
1686      {
1687        return m_shared;
1688    /*  if (m_shared) return true;
1689        if (m_data->isShared())        
1690        {                  
1691            updateShareStatus(true);
1692            return true;
1693        }
1694        return false;*/
1695      }
1696    
1697      void forceResolve()
1698      {
1699        if (isLazy())
1700        {
1701            #ifdef _OPENMP
1702            if (omp_in_parallel())
1703            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1704            throw DataException("Please do not call forceResolve() in a parallel region.");
1705            }
1706            #endif
1707            resolve();
1708        }
1709      }
1710    
1711      /**
1712      \brief if another object is sharing out member data make a copy to work with instead.
1713      This code should only be called from single threaded sections of code.
1714      */
1715      void exclusiveWrite()
1716      {
1717    #ifdef _OPENMP
1718      if (omp_in_parallel())
1719      {
1720    // *((int*)0)=17;
1721        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1722      }
1723    #endif
1724        forceResolve();
1725        if (isShared())
1726        {
1727            DataAbstract* t=m_data->deepCopy();
1728            set_m_data(DataAbstract_ptr(t));
1729        }
1730      }
1731    
1732      /**
1733      \brief checks if caller can have exclusive write to the object
1734      */
1735      void checkExclusiveWrite()
1736      {
1737        if  (isLazy() || isShared())
1738        {
1739            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1740        }
1741      }
1742    
1743      /**
1744      \brief Modify the data abstract hosted by this Data object
1745      For internal use only.
1746      Passing a pointer to null is permitted (do this in the destructor)
1747      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1748      */
1749      void set_m_data(DataAbstract_ptr p);
1750    
1751      friend class DataAbstract;        // To allow calls to updateShareStatus
1752    
1753  };  };
1754    
1755  template <class IValueType>  }   // end namespace escript
1756  void  
1757  Data::initialise(const IValueType& value,  
1758                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1759                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1760    // so that I can dynamic cast between them below.
1761    #include "DataReady.h"
1762    #include "DataLazy.h"
1763    
1764    namespace escript
1765  {  {
1766    //  
1767    // Construct a Data object of the appropriate type.  inline
1768    // Construct the object first as there seems to be a bug which causes  const DataReady*
1769    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1770    // within the shared_ptr constructor.  {
1771    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1772      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1773      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
     m_data=temp_data;  
   } else {  
     DataAbstract* temp=new DataConstant(value,what);  
     boost::shared_ptr<DataAbstract> temp_data(temp);  
     m_data=temp_data;  
   }  
1774  }  }
1775    
1776    inline
1777    DataReady*
1778    Data::getReady()
1779    {
1780       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1781       EsysAssert((dr!=0), "Error - casting to DataReady.");
1782       return dr;
1783    }
1784    
1785    // Be wary of using this for local operations since it (temporarily) increases reference count.
1786    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1787    // getReady() instead
1788    inline
1789    DataReady_ptr
1790    Data::getReadyPtr()
1791    {
1792       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1793       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1794       return dr;
1795    }
1796    
1797    
1798    inline
1799    const_DataReady_ptr
1800    Data::getReadyPtr() const
1801    {
1802       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1803       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1804       return dr;
1805    }
1806    
1807    inline
1808    DataAbstract::ValueType::value_type*
1809    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1810    {
1811       if (isLazy())
1812       {
1813        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1814       }
1815       return getReady()->getSampleDataRW(sampleNo);
1816    }
1817    
1818    inline
1819    const DataAbstract::ValueType::value_type*
1820    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1821    {
1822       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1823       if (l!=0)
1824       {
1825        size_t offset=0;
1826        if (bufferg==NULL)
1827        {
1828            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1829        }
1830        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1831        return &((*res)[offset]);
1832       }
1833       return getReady()->getSampleDataRO(sampleNo);
1834    }
1835    
1836    
1837    
1838    /**
1839       Modify a filename for MPI parallel output to multiple files
1840    */
1841    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1842    
1843  /**  /**
1844     Binary Data object operators.     Binary Data object operators.
1845  */  */
1846  inline double rpow(double x,double y)  inline double rpow(double x,double y)
1847  {  {
1848      return pow(y,x);      return pow(y,x);
1849  };  }
1850    
1851  /**  /**
1852    \brief    \brief
# Line 1514  ESCRIPT_DLL_API Data operator*(const boo Line 1940  ESCRIPT_DLL_API Data operator*(const boo
1940  */  */
1941  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);
1942    
1943    
1944    
1945  /**  /**
1946    \brief    \brief
1947    Output operator    Output operator
# Line 1523  ESCRIPT_DLL_API std::ostream& operator<< Line 1951  ESCRIPT_DLL_API std::ostream& operator<<
1951  /**  /**
1952    \brief    \brief
1953    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1954    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1955    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1956    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1957    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1958  */  */
1959  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1960  Data  Data
1961  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1962                       Data& arg1,                       Data& arg_1,
1963                       int axis_offset=0,                       int axis_offset=0,
1964                       int transpose=0);                       int transpose=0);
1965    
1966  /**  /**
1967    \brief    \brief
   Return true if operands are equivalent, else return false.  
   NB: this operator does very little at this point, and isn't to  
   be relied on. Requires further implementation.  
 */  
 //ESCRIPT_DLL_API bool operator==(const Data& left, const Data& right);  
   
 /**  
   \brief  
1968    Perform the given binary operation with this and right as operands.    Perform the given binary operation with this and right as operands.
1969    Right is a Data object.    Right is a Data object.
1970  */  */
# Line 1556  Data::binaryOp(const Data& right, Line 1976  Data::binaryOp(const Data& right,
1976  {  {
1977     //     //
1978     // 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
1979     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1980       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1981     }     }
1982    
1983       if (isLazy() || right.isLazy())
1984       {
1985         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1986       }
1987     //     //
1988     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1989     Data tempRight(right);     Data tempRight(right);
1990    
1991     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1992       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1993         //         //
1994         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1995         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1996       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1997         //         //
1998         // interpolate onto the RHS function space         // interpolate onto the RHS function space
1999         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2000         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2001           set_m_data(tempLeft.m_data);
2002       }       }
2003     }     }
2004     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1587  Data::binaryOp(const Data& right, Line 2014  Data::binaryOp(const Data& right,
2014       // of any data type       // of any data type
2015       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2016       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2017       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2018     } else if (isTagged()) {     } else if (isTagged()) {
2019       //       //
2020       // Tagged data is operated on serially, the right hand side can be       // Tagged data is operated on serially, the right hand side can be
# Line 1609  Data::binaryOp(const Data& right, Line 2036  Data::binaryOp(const Data& right,
2036       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2037       escript::binaryOp(*leftC,*rightC,operation);       escript::binaryOp(*leftC,*rightC,operation);
2038     }     }
    #if defined DOPROF  
    profData->binary++;  
    #endif  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on other and return the result.  
   Given operation is performed on each element of each data point, thus  
   argument object is a rank n Data object, and returned object is a rank n  
   Data object.  
   Calls Data::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 Data  
 unaryOp(const Data& other,  
         UnaryFunction operation)  
 {  
   Data result;  
   result.copy(other);  
   result.unaryOp(operation);  
   return result;  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on this.  
   Given operation is performed on each element of each data point.  
   Calls escript::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 void  
 Data::unaryOp(UnaryFunction operation)  
 {  
   if (isExpanded()) {  
     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isTagged()) {  
     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isConstant()) {  
     DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");  
     escript::unaryOp(*leftC,operation);  
   }  
2039  }  }
2040    
2041  /**  /**
# Line 1684  Data::algorithm(BinaryFunction operation Line 2062  Data::algorithm(BinaryFunction operation
2062      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2063      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2064      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2065      } else if (isEmpty()) {
2066        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2067      } else if (isLazy()) {
2068        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2069      } else {
2070        throw DataException("Error - Data encapsulates an unknown type.");
2071    }    }
   return 0;  
2072  }  }
2073    
2074  /**  /**
2075    \brief    \brief
2076    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.
2077    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2078    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
2079    rank 0 Data object.    rank 0 Data object.
2080    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2081  */  */
# Line 1701  inline Line 2084  inline
2084  Data  Data
2085  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2086  {  {
2087    if (isExpanded()) {    if (isEmpty()) {
2088      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2089      }
2090      else if (isExpanded()) {
2091        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2092      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2093      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2094      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2095      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2096      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2097      return result;      return result;
2098    } else if (isTagged()) {    }
2099      else if (isTagged()) {
2100      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());  
2101      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2102      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2103        defval[0]=0;
2104        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2105      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2106      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2107    } else if (isConstant()) {    }
2108      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2109        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2110      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2111      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2112      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2113      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2114      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2115      return result;      return result;
2116      } else if (isLazy()) {
2117        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2118      } else {
2119        throw DataException("Error - Data encapsulates an unknown type.");
2120    }    }
2121    Data falseRetVal; // to keep compiler quiet  }
2122    return falseRetVal;  
2123    /**
2124      \brief
2125      Compute a tensor operation with two Data objects
2126      \param arg_0 - Input - Data object
2127      \param arg_1 - Input - Data object
2128      \param operation - Input - Binary op functor
2129    */
2130    template <typename BinaryFunction>
2131    inline
2132    Data
2133    C_TensorBinaryOperation(Data const &arg_0,
2134                            Data const &arg_1,
2135                            BinaryFunction operation)
2136    {
2137      if (arg_0.isEmpty() || arg_1.isEmpty())
2138      {
2139         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2140      }
2141      if (arg_0.isLazy() || arg_1.isLazy())
2142      {
2143         throw DataException("Error - Operations not permitted on lazy data.");
2144      }
2145      // Interpolate if necessary and find an appropriate function space
2146      Data arg_0_Z, arg_1_Z;
2147      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2148        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2149          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2150          arg_1_Z = Data(arg_1);
2151        }
2152        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2153          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2154          arg_0_Z =Data(arg_0);
2155        }
2156        else {
2157          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2158        }
2159      } else {
2160          arg_0_Z = Data(arg_0);
2161          arg_1_Z = Data(arg_1);
2162      }
2163      // Get rank and shape of inputs
2164      int rank0 = arg_0_Z.getDataPointRank();
2165      int rank1 = arg_1_Z.getDataPointRank();
2166      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2167      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2168      int size0 = arg_0_Z.getDataPointSize();
2169      int size1 = arg_1_Z.getDataPointSize();
2170      // Declare output Data object
2171      Data res;
2172    
2173      if (shape0 == shape1) {
2174        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2175          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2176          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2177          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2178          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2179    
2180          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2181        }
2182        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2183    
2184          // Prepare the DataConstant input
2185          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2186    
2187          // Borrow DataTagged input from Data object
2188          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2189    
2190          // Prepare a DataTagged output 2
2191          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2192          res.tag();
2193          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2194    
2195          // Prepare offset into DataConstant
2196          int offset_0 = tmp_0->getPointOffset(0,0);
2197          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2198    
2199          // Get the pointers to the actual data
2200          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2201          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2202    
2203          // Compute a result for the default
2204          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2205          // Compute a result for each tag
2206          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2207          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2208          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2209            tmp_2->addTag(i->first);
2210            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2211            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2212    
2213            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2214          }
2215    
2216        }
2217        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2218          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2219          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2220          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2221          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2222    
2223          int sampleNo_1,dataPointNo_1;
2224          int numSamples_1 = arg_1_Z.getNumSamples();
2225          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2226          int offset_0 = tmp_0->getPointOffset(0,0);
2227          res.requireWrite();
2228          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2229          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2230            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2231              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2232              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2233              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2234              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2235              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2236              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2237            }
2238          }
2239    
2240        }
2241        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2242          // Borrow DataTagged input from Data object
2243          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2244    
2245          // Prepare the DataConstant input
2246          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2247    
2248          // Prepare a DataTagged output 2
2249          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2250          res.tag();
2251          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2252    
2253          // Prepare offset into DataConstant
2254          int offset_1 = tmp_1->getPointOffset(0,0);
2255    
2256          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2257          // Get the pointers to the actual data
2258          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2259          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2260          // Compute a result for the default
2261          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2262          // Compute a result for each tag
2263          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2264          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2265          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2266            tmp_2->addTag(i->first);
2267            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2268            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2269            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2270          }
2271    
2272        }
2273        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2274          // Borrow DataTagged input from Data object
2275          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2276    
2277          // Borrow DataTagged input from Data object
2278          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2279    
2280          // Prepare a DataTagged output 2
2281          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2282          res.tag();        // DataTagged output
2283          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2284    
2285          // Get the pointers to the actual data
2286          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2287          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2288          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2289    
2290          // Compute a result for the default
2291          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2292          // Merge the tags
2293          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2294          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2295          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2296          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2297            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2298          }
2299          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2300            tmp_2->addTag(i->first);
2301          }
2302          // Compute a result for each tag
2303          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2304          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2305    
2306            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2307            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2308            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2309    
2310            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2311          }
2312    
2313        }
2314        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2315          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2316          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2317          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2318          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2319          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2320    
2321          int sampleNo_0,dataPointNo_0;
2322          int numSamples_0 = arg_0_Z.getNumSamples();
2323          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2324          res.requireWrite();
2325          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2326          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2327            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2328            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2329            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2330              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2331              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2332              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2333              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2334              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2335            }
2336          }
2337    
2338        }
2339        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2340          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2341          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2342          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2343          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2344    
2345          int sampleNo_0,dataPointNo_0;
2346          int numSamples_0 = arg_0_Z.getNumSamples();
2347          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2348          int offset_1 = tmp_1->getPointOffset(0,0);
2349          res.requireWrite();
2350          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2351          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2352            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2353              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2354              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2355    
2356              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2357              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2358              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2359    
2360    
2361              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2362            }
2363          }
2364    
2365        }
2366        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2367          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2368          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2369          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2370          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2371          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2372    
2373          int sampleNo_0,dataPointNo_0;
2374          int numSamples_0 = arg_0_Z.getNumSamples();
2375          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2376          res.requireWrite();
2377          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2378          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2379            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2380            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2381            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2382              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2383              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2384              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2385              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2386              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2387            }
2388          }
2389    
2390        }
2391        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2392          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2393          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2394          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2395          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2396          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2397    
2398          int sampleNo_0,dataPointNo_0;
2399          int numSamples_0 = arg_0_Z.getNumSamples();
2400          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2401          res.requireWrite();
2402          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2403          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2404            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2405              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2406              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2407              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2408              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2409              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2410              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2411              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2412            }
2413          }
2414    
2415        }
2416        else {
2417          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2418        }
2419    
2420      } else if (0 == rank0) {
2421        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2422          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2423          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2424          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2425          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2426          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2427        }
2428        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2429    
2430          // Prepare the DataConstant input
2431          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2432    
2433          // Borrow DataTagged input from Data object
2434          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2435    
2436          // Prepare a DataTagged output 2
2437          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2438          res.tag();
2439          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2440    
2441          // Prepare offset into DataConstant
2442          int offset_0 = tmp_0->getPointOffset(0,0);
2443          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2444    
2445          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2446          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2447    
2448          // Compute a result for the default
2449          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2450          // Compute a result for each tag
2451          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2452          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2453          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2454            tmp_2->addTag(i->first);
2455            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2456            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2457            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2458          }
2459    
2460        }
2461        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2462    
2463          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2464          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2465          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2466          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2467    
2468          int sampleNo_1,dataPointNo_1;
2469          int numSamples_1 = arg_1_Z.getNumSamples();
2470          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2471          int offset_0 = tmp_0->getPointOffset(0,0);
2472          res.requireWrite();
2473          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2474          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2475            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2476              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2477              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2478              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2479              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2480              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2481              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2482    
2483            }
2484          }
2485    
2486        }
2487        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2488    
2489          // Borrow DataTagged input from Data object
2490          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2491    
2492          // Prepare the DataConstant input
2493          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2494    
2495          // Prepare a DataTagged output 2
2496          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2497          res.tag();
2498          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2499    
2500          // Prepare offset into DataConstant
2501          int offset_1 = tmp_1->getPointOffset(0,0);
2502          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2503    
2504          // Get the pointers to the actual data
2505          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2506          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2507    
2508    
2509          // Compute a result for the default
2510          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2511          // Compute a result for each tag
2512          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2513          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2514          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2515            tmp_2->addTag(i->first);
2516            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2517            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2518    
2519            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2520          }
2521    
2522        }
2523        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2524    
2525          // Borrow DataTagged input from Data object
2526          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2527    
2528          // Borrow DataTagged input from Data object
2529          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2530    
2531          // Prepare a DataTagged output 2
2532          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2533          res.tag();        // DataTagged output
2534          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2535    
2536          // Get the pointers to the actual data
2537          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2538          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2539          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2540    
2541          // Compute a result for the default
2542          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2543          // Merge the tags
2544          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2545          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2546          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2547          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2548            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2549          }
2550          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2551            tmp_2->addTag(i->first);
2552          }
2553          // Compute a result for each tag
2554          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2555          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2556            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2557            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2558            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2559    
2560            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2561          }
2562    
2563        }
2564        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2565    
2566          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2567          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2568          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2569          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2570          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2571    
2572          int sampleNo_0,dataPointNo_0;
2573          int numSamples_0 = arg_0_Z.getNumSamples();
2574          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2575          res.requireWrite();
2576          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2577          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2578            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2579            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2580            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2581              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2582              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2583              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2584              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2585              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2586            }
2587          }
2588    
2589        }
2590        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2591          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2592          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2593          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2594          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2595    
2596          int sampleNo_0,dataPointNo_0;
2597          int numSamples_0 = arg_0_Z.getNumSamples();
2598          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2599          int offset_1 = tmp_1->getPointOffset(0,0);
2600          res.requireWrite();
2601          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2602          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2603            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2604              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2605              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2606              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2607              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2608              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2609              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2610            }
2611          }
2612    
2613    
2614        }
2615        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2616    
2617          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2618          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2619          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2620          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2621          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2622    
2623          int sampleNo_0,dataPointNo_0;
2624          int numSamples_0 = arg_0_Z.getNumSamples();
2625          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2626          res.requireWrite();
2627          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2628          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2629            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2630            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2631            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2632              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2633              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2634              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2635              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2636              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2637            }
2638          }
2639    
2640        }
2641        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2642    
2643          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2644          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2645          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2646          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2647          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2648    
2649          int sampleNo_0,dataPointNo_0;
2650          int numSamples_0 = arg_0_Z.getNumSamples();
2651          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2652          res.requireWrite();
2653          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2654          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2655            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2656              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2657              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2658              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2659              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2660              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2661              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2662              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2663            }
2664          }
2665    
2666        }
2667        else {
2668          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2669        }
2670    
2671      } else if (0 == rank1) {
2672        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2673          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2674          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2675          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2676          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2677          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2678        }
2679        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2680    
2681          // Prepare the DataConstant input
2682          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2683    
2684          // Borrow DataTagged input from Data object
2685          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2686    
2687          // Prepare a DataTagged output 2
2688          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2689          res.tag();
2690          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2691    
2692          // Prepare offset into DataConstant
2693          int offset_0 = tmp_0->getPointOffset(0,0);
2694          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2695    
2696          //Get the pointers to the actual data
2697          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2698          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2699    
2700          // Compute a result for the default
2701          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2702          // Compute a result for each tag
2703          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2704          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2705          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2706            tmp_2->addTag(i->first);
2707            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2708            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2709            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2710          }
2711        }
2712        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2713    
2714          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2715          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2716          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2717          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2718    
2719          int sampleNo_1,dataPointNo_1;
2720          int numSamples_1 = arg_1_Z.getNumSamples();
2721          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2722          int offset_0 = tmp_0->getPointOffset(0,0);
2723          res.requireWrite();
2724          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2725          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2726            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2727              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2728              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2729              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2730              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2731              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2732              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2733            }
2734          }
2735    
2736        }
2737        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2738    
2739          // Borrow DataTagged input from Data object
2740          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2741    
2742          // Prepare the DataConstant input
2743          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2744    
2745          // Prepare a DataTagged output 2
2746          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2747          res.tag();
2748          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2749    
2750          // Prepare offset into DataConstant
2751          int offset_1 = tmp_1->getPointOffset(0,0);
2752          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2753          // Get the pointers to the actual data
2754          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2755          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2756          // Compute a result for the default
2757          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2758          // Compute a result for each tag
2759          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2760          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2761          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2762            tmp_2->addTag(i->first);
2763            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2764            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2765            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2766          }
2767    
2768        }
2769        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2770    
2771          // Borrow DataTagged input from Data object
2772          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2773    
2774          // Borrow DataTagged input from Data object
2775          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2776    
2777          // Prepare a DataTagged output 2
2778          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2779          res.tag();        // DataTagged output
2780          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2781    
2782          // Get the pointers to the actual data
2783          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2784          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2785          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2786    
2787          // Compute a result for the default
2788          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2789          // Merge the tags
2790          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2791          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2792          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2793          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2794            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2795          }
2796          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2797            tmp_2->addTag(i->first);
2798          }
2799          // Compute a result for each tag
2800          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2801          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
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    
2920      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      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2941      int size0 = arg_0_Z.getDataPointSize();
2942    
2943      // Declare output Data object
2944      Data res;
2945    
2946      if (arg_0_Z.isConstant()) {
2947        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2948        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2949        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2950        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2951      }
2952      else if (arg_0_Z.isTagged()) {
2953    
2954        // Borrow DataTagged input from Data object
2955        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2956    
2957        // Prepare a DataTagged output 2
2958        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2959        res.tag();
2960        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2961    
2962        // Get the pointers to the actual data
2963        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2964        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2965        // Compute a result for the default
2966        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2967        // Compute a result for each tag
2968        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2969        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2970        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2971          tmp_2->addTag(i->first);
2972          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2973          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2974          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2975        }
2976    
2977      }
2978      else if (arg_0_Z.isExpanded()) {
2979    
2980        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2981        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2982        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2983    
2984        int sampleNo_0,dataPointNo_0;
2985        int numSamples_0 = arg_0_Z.getNumSamples();
2986        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2987        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2988        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2989          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2990            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2991            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2992            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2993            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2994            tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2995          }
2996        }
2997      }
2998      else {
2999        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3000      }
3001    
3002      return res;
3003  }  }
3004    
3005  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26