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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26