/[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 2742 by jfenwick, Thu Nov 12 06:03:37 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    
287    ESCRIPT_DLL_API    ESCRIPT_DLL_API
288    const boost::python::numeric::array    void
289    convertToNumArrayFromDPNo(int ProcNo,    setValueOfDataPoint(int dataPointNo, const double);
                                                         int sampleNo,  
                             int dataPointNo);  
 #else  
   ESCRIPT_DLL_API  
   const boost::python::numeric::array  
   convertToNumArrayFromDPNo(int procNo,  
                 int sampleNo,  
                 int dataPointNo);  
 #endif  
   
290    
291    /**    /**
292       \brief       \brief Return a data point across all processors as a python tuple.
      Fills the expanded Data object from values of a python numarray object.  
293    */    */
294    ESCRIPT_DLL_API    ESCRIPT_DLL_API
295    void    const boost::python::object
296    fillFromNumArray(const boost::python::numeric::array);    getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
297    
298    /**    /**
299       \brief       \brief
300       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
301    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
302    */    */
303    ESCRIPT_DLL_API    ESCRIPT_DLL_API
304    int    int
# Line 316  class Data { Line 312  class Data {
312    escriptDataC    escriptDataC
313    getDataC();    getDataC();
314    
315    
316    
317    /**    /**
318       \brief       \brief
319       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 325  class Data { Line 323  class Data {
323    getDataC() const;    getDataC() const;
324    
325    /**    /**
326       \brief      \brief How much space is required to evaulate a sample of the Data.
      Write the data as a string.  
327    */    */
328    ESCRIPT_DLL_API    ESCRIPT_DLL_API
329    inline    size_t
330    std::string    getSampleBufferSize() const;
331    toString() const  
332    {  
     return m_data->toString();  
   }  
333    
334    /**    /**
335       \brief       \brief
336       Return the DataArrayView of the point data. This essentially contains       Write the data as a string. For large amounts of data, a summary is printed.
      the shape information for each data point although it also may be used  
      to manipulate the point data.  
337    */    */
338    ESCRIPT_DLL_API    ESCRIPT_DLL_API
339    inline    std::string
340    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
341    
342    /**    /**
343       \brief       \brief
# Line 363  class Data { Line 352  class Data {
352       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
353       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
354       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
355    */    */
356    ESCRIPT_DLL_API    ESCRIPT_DLL_API
357    void    void
358    tag();    tag();
359    
360    /**    /**
361        \brief If this data is lazy, then convert it to ready data.
362        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
363      */
364      ESCRIPT_DLL_API
365      void
366      resolve();
367    
368    
369      /**
370       \brief Ensures data is ready for write access.
371      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
372      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
373      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
374      Doing so might introduce additional sharing.
375      */
376      ESCRIPT_DLL_API
377      void
378      requireWrite();
379    
380      /**
381       \brief       \brief
382       Return true if this Data is expanded.       Return true if this Data is expanded.
383         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
384    */    */
385    ESCRIPT_DLL_API    ESCRIPT_DLL_API
386    bool    bool
# Line 379  class Data { Line 388  class Data {
388    
389    /**    /**
390       \brief       \brief
391         Return true if this Data is expanded or resolves to expanded.
392         That is, if it has a separate value for each datapoint in the sample.
393      */
394      ESCRIPT_DLL_API
395      bool
396      actsExpanded() const;
397      
398    
399      /**
400         \brief
401       Return true if this Data is tagged.       Return true if this Data is tagged.
402    */    */
403    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 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
819      boost::python::object
820      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    ESCRIPT_DLL_API
828    boost::python::numeric::array    boost::python::object
829    integrate() const;    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 1145  class Data { Line 1354  class Data {
1354    Data& operator/=(const boost::python::object& right);    Data& operator/=(const boost::python::object& right);
1355    
1356    /**    /**
1357        \brief return inverse of matricies.
1358      */
1359      ESCRIPT_DLL_API
1360      Data
1361      matrixInverse() const;
1362    
1363      /**
1364       \brief       \brief
1365       Returns true if this can be interpolated to functionspace.       Returns true if this can be interpolated to functionspace.
1366    */    */
# Line 1203  class Data { Line 1419  class Data {
1419    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1420    inline    inline
1421    void    void
1422    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1423    
1424    /**    /**
1425       \brief       \brief
# Line 1214  class Data { Line 1430  class Data {
1430    */    */
1431    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1432    Data    Data
1433    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1434    
1435    /**    /**
1436       \brief       \brief
# Line 1227  class Data { Line 1443  class Data {
1443    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1444    void    void
1445    setSlice(const Data& value,    setSlice(const Data& value,
1446             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1447    
1448    /**    /**
1449       \brief       \brief
1450       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.  
1451    */    */
1452    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1453    void    void
1454    archiveData(const std::string fileName);          print(void);
1455    
1456    /**    /**
1457       \brief       \brief
1458       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1459       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1460       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.  
1461    */    */
1462    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1463    void          int
1464    extractData(const std::string fileName,          get_MPIRank(void) const;
               const FunctionSpace& fspace);  
   
1465    
1466    /**    /**
1467       \brief       \brief
1468       print the data values to stdout. Used for debugging       return the MPI rank number of the local data
1469                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1470                     is returned
1471    */    */
1472    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1473    void          int
1474      print(void);          get_MPISize(void) const;
1475    
1476    /**    /**
1477       \brief       \brief
1478       return the MPI rank number of the local data       return the MPI rank number of the local data
1479           MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()                   MPI_COMM_WORLD is assumed and returned.
          is returned  
1480    */    */
1481    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1482      int          MPI_Comm
1483      get_MPIRank(void) const;          get_MPIComm(void) const;
1484    
1485    /**    /**
1486       \brief       \brief
1487       return the MPI rank number of the local data       return the object produced by the factory, which is a DataConstant or DataExpanded
1488           MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()      TODO Ownership of this object should be explained in doco.
          is returned  
1489    */    */
1490    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1491      int          DataAbstract*
1492      get_MPISize(void) const;          borrowData(void) const;
1493    
   /**  
      \brief  
      return the MPI rank number of the local data  
          MPI_COMM_WORLD is assumed and returned.  
   */  
1494    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1495      MPI_Comm          DataAbstract_ptr
1496      get_MPIComm(void) const;          borrowDataPtr(void) const;
1497    
1498      ESCRIPT_DLL_API
1499            DataReady_ptr
1500            borrowReadyPtr(void) const;
1501    
1502    
1503    
1504    /**    /**
1505       \brief       \brief
1506       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.
1507         TODO Eventually these should be inlined.
1508         \param i - position(offset) in the underlying datastructure
1509    */    */
1510    
1511      ESCRIPT_DLL_API
1512            DataTypes::ValueType::const_reference
1513            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1514    
1515    
1516      ESCRIPT_DLL_API
1517            DataTypes::ValueType::reference
1518            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1519    
1520    
1521    
1522    /**
1523       \brief Create a buffer for use by getSample
1524       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1525       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1526      
1527       In multi-threaded sections, this needs to be called on each thread.
1528    
1529       \return A BufferGroup* if Data is lazy, NULL otherwise.
1530       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1531    */
1532    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1533      DataAbstract*    BufferGroup*
1534      borrowData(void) const;    allocSampleBuffer() const;
1535    
1536    /**
1537       \brief Free a buffer allocated with allocSampleBuffer.
1538       \param buffer Input - pointer to the buffer to deallocate.
1539    */
1540    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1541    
1542   protected:   protected:
1543    
1544   private:   private:
1545    
1546    template <class BinaryOp>
1547      double
1548    #ifdef PASO_MPI
1549      lazyAlgWorker(double init, MPI_Op mpiop_type);
1550    #else
1551      lazyAlgWorker(double init);
1552    #endif
1553    
1554      double
1555      LsupWorker() const;
1556    
1557      double
1558      supWorker() const;
1559    
1560      double
1561      infWorker() const;
1562    
1563      boost::python::object
1564      integrateWorker() const;
1565    
1566      void
1567      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1568    
1569      void
1570      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1571    
1572      // For internal use in Data.cpp only!
1573      // other uses should call the main entry points and allow laziness
1574      Data
1575      minval_nonlazy() const;
1576    
1577      // For internal use in Data.cpp only!
1578      Data
1579      maxval_nonlazy() const;
1580    
1581    
1582    /**    /**
1583       \brief       \brief
1584       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1372  class Data { Line 1650  class Data {
1650       \brief       \brief
1651       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1652    */    */
1653    template <class IValueType>  
1654    void    void
1655    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1656             const DataTypes::ShapeType& shape,
1657               const FunctionSpace& what,               const FunctionSpace& what,
1658               bool expanded);               bool expanded);
1659    
1660      void
1661      initialise(const WrappedArray& value,
1662                     const FunctionSpace& what,
1663                     bool expanded);
1664    
1665    //    //
1666    // flag to protect the data object against any update    // flag to protect the data object against any update
1667    bool m_protected;    bool m_protected;
1668      mutable bool m_shared;
1669      bool m_lazy;
1670    
1671    //    //
1672    // pointer to the actual data object    // pointer to the actual data object
1673    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1674      DataAbstract_ptr m_data;
1675    
1676    // If possible please use getReadyPtr instead.
1677    // But see warning below.
1678      const DataReady*
1679      getReady() const;
1680    
1681      DataReady*
1682      getReady();
1683    
1684    
1685    // Be wary of using this for local operations since it (temporarily) increases reference count.
1686    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1687    // getReady() instead
1688      DataReady_ptr
1689      getReadyPtr();
1690    
1691      const_DataReady_ptr
1692      getReadyPtr() const;
1693    
1694    
1695      /**
1696       \brief Update the Data's shared flag
1697       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1698       For internal use only.
1699      */
1700      void updateShareStatus(bool nowshared) const
1701      {
1702        m_shared=nowshared;     // m_shared is mutable
1703      }
1704    
1705      // In the isShared() method below:
1706      // A problem would occur if m_data (the address pointed to) were being modified
1707      // while the call m_data->is_shared is being executed.
1708      //
1709      // Q: So why do I think this code can be thread safe/correct?
1710      // A: We need to make some assumptions.
1711      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1712      // 2. We assume that no constructions or assignments which will share previously unshared
1713      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1714    //    //
1715    // pointer to the internal profiling data    // This means that the only transition we need to consider, is when a previously shared object is
1716    struct profDataEntry *profData;    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1717      // In those cases the m_shared flag changes to false after m_data has completed changing.
1718      // For any threads executing before the flag switches they will assume the object is still shared.
1719      bool isShared() const
1720      {
1721        return m_shared;
1722    /*  if (m_shared) return true;
1723        if (m_data->isShared())        
1724        {                  
1725            updateShareStatus(true);
1726            return true;
1727        }
1728        return false;*/
1729      }
1730    
1731      void forceResolve()
1732      {
1733        if (isLazy())
1734        {
1735            #ifdef _OPENMP
1736            if (omp_in_parallel())
1737            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1738            throw DataException("Please do not call forceResolve() in a parallel region.");
1739            }
1740            #endif
1741            resolve();
1742        }
1743      }
1744    
1745      /**
1746      \brief if another object is sharing out member data make a copy to work with instead.
1747      This code should only be called from single threaded sections of code.
1748      */
1749      void exclusiveWrite()
1750      {
1751    #ifdef _OPENMP
1752      if (omp_in_parallel())
1753      {
1754    // *((int*)0)=17;
1755        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1756      }
1757    #endif
1758        forceResolve();
1759        if (isShared())
1760        {
1761            DataAbstract* t=m_data->deepCopy();
1762            set_m_data(DataAbstract_ptr(t));
1763        }
1764      }
1765    
1766      /**
1767      \brief checks if caller can have exclusive write to the object
1768      */
1769      void checkExclusiveWrite()
1770      {
1771        if  (isLazy() || isShared())
1772        {
1773            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1774        }
1775      }
1776    
1777      /**
1778      \brief Modify the data abstract hosted by this Data object
1779      For internal use only.
1780      Passing a pointer to null is permitted (do this in the destructor)
1781      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1782      */
1783      void set_m_data(DataAbstract_ptr p);
1784    
1785      friend class DataAbstract;        // To allow calls to updateShareStatus
1786    
1787  };  };
1788    
1789  template <class IValueType>  }   // end namespace escript
1790  void  
1791  Data::initialise(const IValueType& value,  
1792                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1793                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1794    // so that I can dynamic cast between them below.
1795    #include "DataReady.h"
1796    #include "DataLazy.h"
1797    
1798    namespace escript
1799  {  {
1800    //  
1801    // Construct a Data object of the appropriate type.  inline
1802    // Construct the object first as there seems to be a bug which causes  const DataReady*
1803    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1804    // within the shared_ptr constructor.  {
1805    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1806      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1807      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
     m_data=temp_data;  
   } else {  
     DataAbstract* temp=new DataConstant(value,what);  
     boost::shared_ptr<DataAbstract> temp_data(temp);  
     m_data=temp_data;  
   }  
1808  }  }
1809    
1810    inline
1811    DataReady*
1812    Data::getReady()
1813    {
1814       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1815       EsysAssert((dr!=0), "Error - casting to DataReady.");
1816       return dr;
1817    }
1818    
1819    // Be wary of using this for local operations since it (temporarily) increases reference count.
1820    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1821    // getReady() instead
1822    inline
1823    DataReady_ptr
1824    Data::getReadyPtr()
1825    {
1826       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1827       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1828       return dr;
1829    }
1830    
1831    
1832    inline
1833    const_DataReady_ptr
1834    Data::getReadyPtr() const
1835    {
1836       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1837       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1838       return dr;
1839    }
1840    
1841    inline
1842    DataAbstract::ValueType::value_type*
1843    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1844    {
1845       if (isLazy())
1846       {
1847        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1848       }
1849       return getReady()->getSampleDataRW(sampleNo);
1850    }
1851    
1852    inline
1853    const DataAbstract::ValueType::value_type*
1854    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1855    {
1856       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1857       if (l!=0)
1858       {
1859        size_t offset=0;
1860        if (bufferg==NULL)
1861        {
1862            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1863        }
1864        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1865        return &((*res)[offset]);
1866       }
1867       return getReady()->getSampleDataRO(sampleNo);
1868    }
1869    
1870    
1871    
1872    /**
1873       Modify a filename for MPI parallel output to multiple files
1874    */
1875    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1876    
1877  /**  /**
1878     Binary Data object operators.     Binary Data object operators.
1879  */  */
1880  inline double rpow(double x,double y)  inline double rpow(double x,double y)
1881  {  {
1882      return pow(y,x);      return pow(y,x);
1883  };  }
1884    
1885  /**  /**
1886    \brief    \brief
# Line 1514  ESCRIPT_DLL_API Data operator*(const boo Line 1974  ESCRIPT_DLL_API Data operator*(const boo
1974  */  */
1975  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);
1976    
1977    
1978    
1979  /**  /**
1980    \brief    \brief
1981    Output operator    Output operator
# Line 1523  ESCRIPT_DLL_API std::ostream& operator<< Line 1985  ESCRIPT_DLL_API std::ostream& operator<<
1985  /**  /**
1986    \brief    \brief
1987    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1988    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1989    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1990    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1991    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1992  */  */
1993  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1994  Data  Data
1995  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1996                       Data& arg1,                       Data& arg_1,
1997                       int axis_offset=0,                       int axis_offset=0,
1998                       int transpose=0);                       int transpose=0);
1999    
2000  /**  /**
2001    \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  
2002    Perform the given binary operation with this and right as operands.    Perform the given binary operation with this and right as operands.
2003    Right is a Data object.    Right is a Data object.
2004  */  */
# Line 1556  Data::binaryOp(const Data& right, Line 2010  Data::binaryOp(const Data& right,
2010  {  {
2011     //     //
2012     // 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
2013     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
2014       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.");
2015     }     }
2016    
2017       if (isLazy() || right.isLazy())
2018       {
2019         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
2020       }
2021     //     //
2022     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
2023     Data tempRight(right);     Data tempRight(right);
2024    
2025     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
2026       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
2027         //         //
2028         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
2029         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
2030       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
2031         //         //
2032         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2033         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2034         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2035           set_m_data(tempLeft.m_data);
2036       }       }
2037     }     }
2038     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1587  Data::binaryOp(const Data& right, Line 2048  Data::binaryOp(const Data& right,
2048       // of any data type       // of any data type
2049       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2050       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2051       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2052     } else if (isTagged()) {     } else if (isTagged()) {
2053       //       //
2054       // 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 2070  Data::binaryOp(const Data& right,
2070       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2071       escript::binaryOp(*leftC,*rightC,operation);       escript::binaryOp(*leftC,*rightC,operation);
2072     }     }
    #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);  
   }  
2073  }  }
2074    
2075  /**  /**
# Line 1684  Data::algorithm(BinaryFunction operation Line 2096  Data::algorithm(BinaryFunction operation
2096      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2097      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2098      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2099      } else if (isEmpty()) {
2100        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2101      } else if (isLazy()) {
2102        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2103      } else {
2104        throw DataException("Error - Data encapsulates an unknown type.");
2105    }    }
   return 0;  
2106  }  }
2107    
2108  /**  /**
2109    \brief    \brief
2110    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.
2111    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2112    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
2113    rank 0 Data object.    rank 0 Data object.
2114    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2115  */  */
# Line 1701  inline Line 2118  inline
2118  Data  Data
2119  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2120  {  {
2121    if (isExpanded()) {    if (isEmpty()) {
2122      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2123      }
2124      else if (isExpanded()) {
2125        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2126      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2127      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2128      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2129      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2130      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2131      return result;      return result;
2132    } else if (isTagged()) {    }
2133      else if (isTagged()) {
2134      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());  
2135      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2136      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2137        defval[0]=0;
2138        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2139      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2140      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2141    } else if (isConstant()) {    }
2142      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2143        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2144      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2145      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2146      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2147      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2148      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2149      return result;      return result;
2150      } else if (isLazy()) {
2151        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2152      } else {
2153        throw DataException("Error - Data encapsulates an unknown type.");
2154      }
2155    }
2156    
2157    /**
2158      \brief
2159      Compute a tensor operation with two Data objects
2160      \param arg_0 - Input - Data object
2161      \param arg_1 - Input - Data object
2162      \param operation - Input - Binary op functor
2163    */
2164    template <typename BinaryFunction>
2165    inline
2166    Data
2167    C_TensorBinaryOperation(Data const &arg_0,
2168                            Data const &arg_1,
2169                            BinaryFunction operation)
2170    {
2171      if (arg_0.isEmpty() || arg_1.isEmpty())
2172      {
2173         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2174      }
2175      if (arg_0.isLazy() || arg_1.isLazy())
2176      {
2177         throw DataException("Error - Operations not permitted on lazy data.");
2178      }
2179      // Interpolate if necessary and find an appropriate function space
2180      Data arg_0_Z, arg_1_Z;
2181      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2182        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2183          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2184          arg_1_Z = Data(arg_1);
2185        }
2186        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2187          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2188          arg_0_Z =Data(arg_0);
2189        }
2190        else {
2191          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2192        }
2193      } else {
2194          arg_0_Z = Data(arg_0);
2195          arg_1_Z = Data(arg_1);
2196      }
2197      // Get rank and shape of inputs
2198      int rank0 = arg_0_Z.getDataPointRank();
2199      int rank1 = arg_1_Z.getDataPointRank();
2200      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2201      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2202      int size0 = arg_0_Z.getDataPointSize();
2203      int size1 = arg_1_Z.getDataPointSize();
2204      // Declare output Data object
2205      Data res;
2206    
2207      if (shape0 == shape1) {
2208        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2209          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2210          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2211          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2212          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2213    
2214          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2215        }
2216        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2217    
2218          // Prepare the DataConstant input
2219          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2220    
2221          // Borrow DataTagged input from Data object
2222          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2223    
2224          // Prepare a DataTagged output 2
2225          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2226          res.tag();
2227          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2228    
2229          // Prepare offset into DataConstant
2230          int offset_0 = tmp_0->getPointOffset(0,0);
2231          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2232    
2233          // Get the pointers to the actual data
2234          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2235          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2236    
2237          // Compute a result for the default
2238          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2239          // Compute a result for each tag
2240          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2241          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2242          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2243            tmp_2->addTag(i->first);
2244            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2245            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2246    
2247            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2248          }
2249    
2250        }
2251        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2252          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2253          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2254          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2255          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2256    
2257          int sampleNo_1,dataPointNo_1;
2258          int numSamples_1 = arg_1_Z.getNumSamples();
2259          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2260          int offset_0 = tmp_0->getPointOffset(0,0);
2261          res.requireWrite();
2262          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2263          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2264            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2265              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2266              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2267              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2268              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2269              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2270              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2271            }
2272          }
2273    
2274        }
2275        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2276          // Borrow DataTagged input from Data object
2277          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2278    
2279          // Prepare the DataConstant input
2280          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2281    
2282          // Prepare a DataTagged output 2
2283          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2284          res.tag();
2285          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2286    
2287          // Prepare offset into DataConstant
2288          int offset_1 = tmp_1->getPointOffset(0,0);
2289    
2290          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2291          // Get the pointers to the actual data
2292          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2293          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2294          // Compute a result for the default
2295          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2296          // Compute a result for each tag
2297          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2298          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2299          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2300            tmp_2->addTag(i->first);
2301            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2302            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2303            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2304          }
2305    
2306        }
2307        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2308          // Borrow DataTagged input from Data object
2309          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2310    
2311          // Borrow DataTagged input from Data object
2312          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2313    
2314          // Prepare a DataTagged output 2
2315          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2316          res.tag();        // DataTagged output
2317          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2318    
2319          // Get the pointers to the actual data
2320          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2321          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2322          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2323    
2324          // Compute a result for the default
2325          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2326          // Merge the tags
2327          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2328          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2329          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2330          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2331            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2332          }
2333          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2334            tmp_2->addTag(i->first);
2335          }
2336          // Compute a result for each tag
2337          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2338          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2339    
2340            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2341            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2342            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2343    
2344            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2345          }
2346    
2347        }
2348        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2349          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2350          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2351          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2352          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(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          res.requireWrite();
2359          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2360          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2361            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2362            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2363            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2364              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2365              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2366              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2367              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2368              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2369            }
2370          }
2371    
2372        }
2373        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2374          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2375          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2376          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2377          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2378    
2379          int sampleNo_0,dataPointNo_0;
2380          int numSamples_0 = arg_0_Z.getNumSamples();
2381          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2382          int offset_1 = tmp_1->getPointOffset(0,0);
2383          res.requireWrite();
2384          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2385          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2386            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2387              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2388              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2389    
2390              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2391              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2392              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2393    
2394    
2395              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2396            }
2397          }
2398    
2399        }
2400        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2401          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2402          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2403          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2404          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2405          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2406    
2407          int sampleNo_0,dataPointNo_0;
2408          int numSamples_0 = arg_0_Z.getNumSamples();
2409          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2410          res.requireWrite();
2411          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2412          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2413            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2414            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2415            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2416              int offset_0 = tmp_0->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              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2420              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2421            }
2422          }
2423    
2424        }
2425        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2426          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2427          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2428          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2429          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2430          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2431    
2432          int sampleNo_0,dataPointNo_0;
2433          int numSamples_0 = arg_0_Z.getNumSamples();
2434          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2435          res.requireWrite();
2436          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2437          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2438          dataPointNo_0=0;
2439    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2440              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2441              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2442              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2443              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2444              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2445              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2446              tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2447    //       }
2448          }
2449    
2450        }
2451        else {
2452          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2453        }
2454    
2455      } else if (0 == rank0) {
2456        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2457          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2458          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2459          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2460          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2461          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2462        }
2463        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2464    
2465          // Prepare the DataConstant input
2466          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2467    
2468          // Borrow DataTagged input from Data object
2469          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2470    
2471          // Prepare a DataTagged output 2
2472          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2473          res.tag();
2474          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2475    
2476          // Prepare offset into DataConstant
2477          int offset_0 = tmp_0->getPointOffset(0,0);
2478          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2479    
2480          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2481          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2482    
2483          // Compute a result for the default
2484          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2485          // Compute a result for each tag
2486          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2487          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2488          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2489            tmp_2->addTag(i->first);
2490            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2491            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2492            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2493          }
2494    
2495        }
2496        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2497    
2498          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2499          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2500          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2501          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2502    
2503          int sampleNo_1,dataPointNo_1;
2504          int numSamples_1 = arg_1_Z.getNumSamples();
2505          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2506          int offset_0 = tmp_0->getPointOffset(0,0);
2507          res.requireWrite();
2508          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2509          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2510            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2511              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2512              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2513              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2514              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2515              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2516              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2517    
2518            }
2519          }
2520    
2521        }
2522        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2523    
2524          // Borrow DataTagged input from Data object
2525          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2526    
2527          // Prepare the DataConstant input
2528          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2529    
2530          // Prepare a DataTagged output 2
2531          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2532          res.tag();
2533          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2534    
2535          // Prepare offset into DataConstant
2536          int offset_1 = tmp_1->getPointOffset(0,0);
2537          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2538    
2539          // Get the pointers to the actual data
2540          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2541          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2542    
2543    
2544          // Compute a result for the default
2545          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2546          // Compute a result for each tag
2547          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2548          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2549          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2550            tmp_2->addTag(i->first);
2551            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2552            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2553    
2554            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2555          }
2556    
2557        }
2558        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2559    
2560          // Borrow DataTagged input from Data object
2561          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2562    
2563          // Borrow DataTagged input from Data object
2564          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2565    
2566          // Prepare a DataTagged output 2
2567          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2568          res.tag();        // DataTagged output
2569          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2570    
2571          // Get the pointers to the actual data
2572          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2573          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2574          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2575    
2576          // Compute a result for the default
2577          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2578          // Merge the tags
2579          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2580          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2581          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2582          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2583            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2584          }
2585          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2586            tmp_2->addTag(i->first);
2587          }
2588          // Compute a result for each tag
2589          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2590          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2591            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2592            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2593            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2594    
2595            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2596          }
2597    
2598        }
2599        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2600    
2601          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2602          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2603          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2604          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2605          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2606    
2607          int sampleNo_0,dataPointNo_0;
2608          int numSamples_0 = arg_0_Z.getNumSamples();
2609          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
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            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2614            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2615            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2616              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2617              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2618              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2619              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2620              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2621            }
2622          }
2623    
2624        }
2625        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2626          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2627          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2628          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2629          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2630    
2631          int sampleNo_0,dataPointNo_0;
2632          int numSamples_0 = arg_0_Z.getNumSamples();
2633          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2634          int offset_1 = tmp_1->getPointOffset(0,0);
2635          res.requireWrite();
2636          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2637          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2638            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2639              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2640              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2641              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2642              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2643              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2644              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2645            }
2646          }
2647    
2648    
2649        }
2650        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2651    
2652          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2653          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2654          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2655          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2656          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2657    
2658          int sampleNo_0,dataPointNo_0;
2659          int numSamples_0 = arg_0_Z.getNumSamples();
2660          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2661          res.requireWrite();
2662          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2663          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2664            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2665            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2666            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2667              int offset_0 = tmp_0->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              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2671              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2672            }
2673          }
2674    
2675        }
2676        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2677    
2678          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2679          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2680          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2681          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2682          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2683    
2684          int sampleNo_0,dataPointNo_0;
2685          int numSamples_0 = arg_0_Z.getNumSamples();
2686          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2687          res.requireWrite();
2688          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2689          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2690            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2691              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2692              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2693              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2694              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2695              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2696              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2697              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2698            }
2699          }
2700    
2701        }
2702        else {
2703          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2704        }
2705    
2706      } else if (0 == rank1) {
2707        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2708          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2709          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2710          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2711          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2712          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2713        }
2714        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2715    
2716          // Prepare the DataConstant input
2717          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2718    
2719          // Borrow DataTagged input from Data object
2720          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2721    
2722          // Prepare a DataTagged output 2
2723          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2724          res.tag();
2725          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2726    
2727          // Prepare offset into DataConstant
2728          int offset_0 = tmp_0->getPointOffset(0,0);
2729          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2730    
2731          //Get the pointers to the actual data
2732          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2733          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2734    
2735          // Compute a result for the default
2736          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2737          // Compute a result for each tag
2738          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2739          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2740          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2741            tmp_2->addTag(i->first);
2742            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2743            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2744            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2745          }
2746        }
2747        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2748    
2749          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2750          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2751          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2752          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2753    
2754          int sampleNo_1,dataPointNo_1;
2755          int numSamples_1 = arg_1_Z.getNumSamples();
2756          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2757          int offset_0 = tmp_0->getPointOffset(0,0);
2758          res.requireWrite();
2759          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2760          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2761            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2762              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2763              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2764              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2765              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2766              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2767              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2768            }
2769          }
2770    
2771        }
2772        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2773    
2774          // Borrow DataTagged input from Data object
2775          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2776    
2777          // Prepare the DataConstant input
2778          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2779    
2780          // Prepare a DataTagged output 2
2781          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2782          res.tag();
2783          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2784    
2785          // Prepare offset into DataConstant
2786          int offset_1 = tmp_1->getPointOffset(0,0);
2787          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2788          // Get the pointers to the actual data
2789          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2790          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2791          // Compute a result for the default
2792          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2793          // Compute a result for each tag
2794          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2795          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2796          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2797            tmp_2->addTag(i->first);
2798            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2799            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2800            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2801          }
2802    
2803        }
2804        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2805    
2806          // Borrow DataTagged input from Data object
2807          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2808    
2809          // Borrow DataTagged input from Data object
2810          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2811    
2812          // Prepare a DataTagged output 2
2813          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2814          res.tag();        // DataTagged output
2815          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2816    
2817          // Get the pointers to the actual data
2818          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2819          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2820          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2821    
2822          // Compute a result for the default
2823          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2824          // Merge the tags
2825          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2826          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2827          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2828          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2829            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2830          }
2831          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2832            tmp_2->addTag(i->first);
2833          }
2834          // Compute a result for each tag
2835          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2836          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2837            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2838            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2839            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2840            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2841          }
2842    
2843        }
2844        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2845    
2846          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2847          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2848          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2849          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2850          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2851    
2852          int sampleNo_0,dataPointNo_0;
2853          int numSamples_0 = arg_0_Z.getNumSamples();
2854          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
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            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2859            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2860            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2861              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2862              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2863              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2864              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2865              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2866            }
2867          }
2868    
2869        }
2870        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2871          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2872          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2873          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2874          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2875    
2876          int sampleNo_0,dataPointNo_0;
2877          int numSamples_0 = arg_0_Z.getNumSamples();
2878          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2879          int offset_1 = tmp_1->getPointOffset(0,0);
2880          res.requireWrite();
2881          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2882          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2883            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2884              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2885              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2886              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2887              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2888              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2889              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2890            }
2891          }
2892    
2893    
2894        }
2895        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2896    
2897          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2898          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2899          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2900          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2901          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2902    
2903          int sampleNo_0,dataPointNo_0;
2904          int numSamples_0 = arg_0_Z.getNumSamples();
2905          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2906          res.requireWrite();
2907          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2908          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2909            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2910            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2911            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2912              int offset_0 = tmp_0->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              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2916              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2917            }
2918          }
2919    
2920        }
2921        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2922    
2923          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2924          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2925          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2926          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2927          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2928    
2929          int sampleNo_0,dataPointNo_0;
2930          int numSamples_0 = arg_0_Z.getNumSamples();
2931          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2932          res.requireWrite();
2933          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2934          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2935            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2936              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2937              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2938              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2939              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2940              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2941              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2942              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2943            }
2944          }
2945    
2946        }
2947        else {
2948          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2949        }
2950    
2951      } else {
2952        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2953    }    }
2954    Data falseRetVal; // to keep compiler quiet  
2955    return falseRetVal;    return res;
2956    }
2957    
2958    template <typename UnaryFunction>
2959    Data
2960    C_TensorUnaryOperation(Data const &arg_0,
2961                           UnaryFunction operation)
2962    {
2963      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2964      {
2965         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2966      }
2967      if (arg_0.isLazy())
2968      {
2969         throw DataException("Error - Operations not permitted on lazy data.");
2970      }
2971      // Interpolate if necessary and find an appropriate function space
2972      Data arg_0_Z = Data(arg_0);
2973    
2974      // Get rank and shape of inputs
2975      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2976      int size0 = arg_0_Z.getDataPointSize();
2977    
2978      // Declare output Data object
2979      Data res;
2980    
2981      if (arg_0_Z.isConstant()) {
2982        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2983        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2984        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2985        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2986      }
2987      else if (arg_0_Z.isTagged()) {
2988    
2989        // Borrow DataTagged input from Data object
2990        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2991    
2992        // Prepare a DataTagged output 2
2993        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2994        res.tag();
2995        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2996    
2997        // Get the pointers to the actual data
2998        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2999        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
3000        // Compute a result for the default
3001        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3002        // Compute a result for each tag
3003        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3004        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3005        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3006          tmp_2->addTag(i->first);
3007          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
3008          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
3009          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3010        }
3011    
3012      }
3013      else if (arg_0_Z.isExpanded()) {
3014    
3015        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
3016        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3017        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3018    
3019        int sampleNo_0,dataPointNo_0;
3020        int numSamples_0 = arg_0_Z.getNumSamples();
3021        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3022        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3023        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3024        dataPointNo_0=0;
3025    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3026            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3027            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3028            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3029            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3030            tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3031    //      }
3032        }
3033      }
3034      else {
3035        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3036      }
3037    
3038      return res;
3039  }  }
3040    
3041  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26