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

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

  ViewVC Help
Powered by ViewVC 1.1.26