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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26