/[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

temp/escript/src/Data.h revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/escript/src/Data.h revision 2721 by jfenwick, Fri Oct 16 05:40:12 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /** \file Data.h */  /** \file Data.h */
16    
# Line 19  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 26  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" doesn't belong in this file...causes trouble for BruceFactory.cpp */  //#include <omp.h>
33  }  }
34    
35  #include "esysmpi.h"  #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>  
44    #include "BufferGroup.h"
45    
46  namespace escript {  namespace escript {
47    
# Line 47  namespace escript { Line 50  namespace escript {
50  class DataConstant;  class DataConstant;
51  class DataTagged;  class DataTagged;
52  class DataExpanded;  class DataExpanded;
53    class DataLazy;
54    
55  /**  /**
56     \brief     \brief
57     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
58    
59     Description:     Description:
60     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
61     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
62     the object created for the lifetime of the object.     of the Data object.
63     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
64     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
65       Doing so will lead to invalid memory access.
66       This should not affect any methods exposed via boost::python.
67  */  */
68  class Data {  class Data {
69    
# Line 101  class Data { Line 106  class Data {
106         const FunctionSpace& what);         const FunctionSpace& what);
107    
108    /**    /**
109       \brief      \brief Copy Data from an existing vector
110       Constructor which copies data from a DataArrayView.    */
111    
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
112    ESCRIPT_DLL_API    ESCRIPT_DLL_API
113    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
114         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
115         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
116                     bool expanded=false);
117    
118    /**    /**
119       \brief       \brief
120       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
121    
122       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
123       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 128  class Data { Line 128  class Data {
128    */    */
129    ESCRIPT_DLL_API    ESCRIPT_DLL_API
130    Data(double value,    Data(double value,
131         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
132         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
133         bool expanded=false);         bool expanded=false);
134    
# Line 141  class Data { Line 141  class Data {
141    */    */
142    ESCRIPT_DLL_API    ESCRIPT_DLL_API
143    Data(const Data& inData,    Data(const Data& inData,
144         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
145    
146    /**    /**
147       \brief       \brief
148       Constructor which will create Tagged data if expanded is false.       Constructor which copies data from any object that can be treated like a python array/sequence.
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
     ==>*  
   */  
   ESCRIPT_DLL_API  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   ESCRIPT_DLL_API  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from any object that can be converted into  
      a python numarray.  
149    
150       \param value - Input - Input data.       \param value - Input - Input data.
151       \param what - Input - A description of what this data represents.       \param what - Input - A description of what this data represents.
# Line 198  class Data { Line 161  class Data {
161    /**    /**
162       \brief       \brief
163       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
164       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
165       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
166    
167       \param value - Input - Input data.       \param value - Input - Input data.
168       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
# Line 217  class Data { Line 180  class Data {
180         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
181         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
182         bool expanded=false);         bool expanded=false);
183    
184    
185    
186      /**
187        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
188        Once you have passed the pointer, do not delete it.
189      */
190      ESCRIPT_DLL_API
191      explicit Data(DataAbstract* underlyingdata);
192    
193      /**
194        \brief Create a Data based on the supplied DataAbstract
195      */
196      ESCRIPT_DLL_API
197      explicit Data(DataAbstract_ptr underlyingdata);
198    
199    /**    /**
200       \brief       \brief
201       Destructor       Destructor
# Line 225  class Data { Line 204  class Data {
204    ~Data();    ~Data();
205    
206    /**    /**
207       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
208    */    */
209    ESCRIPT_DLL_API    ESCRIPT_DLL_API
210    void    void
211    copy(const Data& other);    copy(const Data& other);
212    
213    /**    /**
214         \brief Return a pointer to a deep copy of this object.
215      */
216      ESCRIPT_DLL_API
217      Data
218      copySelf();
219    
220    
221      /**
222         \brief produce a delayed evaluation version of this Data.
223      */
224      ESCRIPT_DLL_API
225      Data
226      delay();
227    
228      /**
229         \brief convert the current data into lazy data.
230      */
231      ESCRIPT_DLL_API
232      void
233      delaySelf();
234    
235    
236      /**
237       Member access methods.       Member access methods.
238    */    */
239    
# Line 247  class Data { Line 248  class Data {
248    
249    /**    /**
250       \brief       \brief
251       Returns trueif the data object is protected against update       Returns true, if the data object is protected against update
252    
253    */    */
254    ESCRIPT_DLL_API    ESCRIPT_DLL_API
255    bool    bool
256    isProtected() const;    isProtected() const;
257    
258    /**  
259       \brief  /**
260       Return the values of a data point on this process     \brief
261    */     Return the value of a data point as a python tuple.
262    */
263    ESCRIPT_DLL_API    ESCRIPT_DLL_API
264    const boost::python::numeric::array    const boost::python::object
265    getValueOfDataPoint(int dataPointNo);    getValueOfDataPointAsTuple(int dataPointNo);
266    
267    /**    /**
268       \brief       \brief
# Line 272  class Data { Line 274  class Data {
274    
275    /**    /**
276       \brief       \brief
277       sets the values of a data-point from a numarray object on this process       sets the values of a data-point from a array-like object on this process
278    */    */
279    ESCRIPT_DLL_API    ESCRIPT_DLL_API
280    void    void
281    setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array&);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282    
283    /**    /**
284       \brief       \brief
# Line 287  class Data { Line 289  class Data {
289    setValueOfDataPoint(int dataPointNo, const double);    setValueOfDataPoint(int dataPointNo, const double);
290    
291    /**    /**
292       \brief       \brief Return a data point across all processors as a python tuple.
      Return the value of the specified data-point across all processors  
293    */    */
294    ESCRIPT_DLL_API    ESCRIPT_DLL_API
295    const boost::python::numeric::array    const boost::python::object
296    getValueOfGlobalDataPoint(int procNo, int dataPointNo);    getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
297    
298    /**    /**
299       \brief       \brief
300       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
301    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
302    */    */
303    ESCRIPT_DLL_API    ESCRIPT_DLL_API
304    int    int
# Line 313  class Data { Line 312  class Data {
312    escriptDataC    escriptDataC
313    getDataC();    getDataC();
314    
315    
316    
317    /**    /**
318       \brief       \brief
319       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 322  class Data { Line 323  class Data {
323    getDataC() const;    getDataC() const;
324    
325    /**    /**
326       \brief      \brief How much space is required to evaulate a sample of the Data.
      Write the data as a string.  
327    */    */
328    ESCRIPT_DLL_API    ESCRIPT_DLL_API
329    inline    size_t
330    std::string    getSampleBufferSize() const;
331    toString() const  
332    {  
     return m_data->toString();  
   }  
333    
334    /**    /**
335       \brief       \brief
336       Return the DataArrayView of the point data. This essentially contains       Write the data as a string. For large amounts of data, a summary is printed.
      the shape information for each data point although it also may be used  
      to manipulate the point data.  
337    */    */
338    ESCRIPT_DLL_API    ESCRIPT_DLL_API
339    inline    std::string
340    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
341    
342    /**    /**
343       \brief       \brief
# Line 360  class Data { Line 352  class Data {
352       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
353       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
354       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
355    */    */
356    ESCRIPT_DLL_API    ESCRIPT_DLL_API
357    void    void
358    tag();    tag();
359    
360    /**    /**
361        \brief If this data is lazy, then convert it to ready data.
362        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
363      */
364      ESCRIPT_DLL_API
365      void
366      resolve();
367    
368    
369      /**
370       \brief Ensures data is ready for write access.
371      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
372      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
373      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
374      Doing so might introduce additional sharing.
375      */
376      ESCRIPT_DLL_API
377      void
378      requireWrite();
379    
380      /**
381       \brief       \brief
382       Return true if this Data is expanded.       Return true if this Data is expanded.
383         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
384    */    */
385    ESCRIPT_DLL_API    ESCRIPT_DLL_API
386    bool    bool
# Line 376  class Data { Line 388  class Data {
388    
389    /**    /**
390       \brief       \brief
391         Return true if this Data is expanded or resolves to expanded.
392         That is, if it has a separate value for each datapoint in the sample.
393      */
394      ESCRIPT_DLL_API
395      bool
396      actsExpanded() const;
397      
398    
399      /**
400         \brief
401       Return true if this Data is tagged.       Return true if this Data is tagged.
402    */    */
403    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 391  class Data { Line 413  class Data {
413    isConstant() const;    isConstant() const;
414    
415    /**    /**
416         \brief Return true if this Data is lazy.
417      */
418      ESCRIPT_DLL_API
419      bool
420      isLazy() const;
421    
422      /**
423         \brief Return true if this data is ready.
424      */
425      ESCRIPT_DLL_API
426      bool
427      isReady() const;
428    
429      /**
430       \brief       \brief
431       Return true if this Data is empty.       Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
432    contains datapoints.
433    */    */
434    ESCRIPT_DLL_API    ESCRIPT_DLL_API
435    bool    bool
# Line 424  class Data { Line 461  class Data {
461    */    */
462    ESCRIPT_DLL_API    ESCRIPT_DLL_API
463    inline    inline
464    const AbstractDomain&  //   const AbstractDomain&
465      const_Domain_ptr
466    getDomain() const    getDomain() const
467    {    {
468       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
469    }    }
470    
471    
472      /**
473         \brief
474         Return the domain.
475         TODO: For internal use only.   This should be removed.
476      */
477      ESCRIPT_DLL_API
478      inline
479    //   const AbstractDomain&
480      Domain_ptr
481      getDomainPython() const
482      {
483         return getFunctionSpace().getDomainPython();
484      }
485    
486    /**    /**
487       \brief       \brief
488       Return a copy of the domain.       Return a copy of the domain.
# Line 444  class Data { Line 497  class Data {
497    */    */
498    ESCRIPT_DLL_API    ESCRIPT_DLL_API
499    inline    inline
500    int    unsigned int
501    getDataPointRank() const    getDataPointRank() const
502    {    {
503      return m_data->getPointDataView().getRank();      return m_data->getRank();
504    }    }
505    
506    /**    /**
# Line 484  class Data { Line 537  class Data {
537    {    {
538      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
539    }    }
540    
541    
542      /**
543        \brief
544        Return the number of values in the shape for this object.
545      */
546      ESCRIPT_DLL_API
547      int
548      getNoValues() const
549      {
550        return m_data->getNoValues();
551      }
552    
553    
554    /**    /**
555       \brief       \brief
556       dumps the object into a netCDF file       dumps the object into a netCDF file
# Line 491  class Data { Line 558  class Data {
558    ESCRIPT_DLL_API    ESCRIPT_DLL_API
559    void    void
560    dump(const std::string fileName) const;    dump(const std::string fileName) const;
561    
562     /**
563      \brief returns the values of the object as a list of tuples (one for each datapoint).
564    
565      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
566    If false, the result is a list of scalars [1, 2, ...]
567     */
568      ESCRIPT_DLL_API
569      const boost::python::object
570      toListOfTuples(bool scalarastuple=true);
571    
572    
573     /**
574        \brief
575        Return the sample data for the given sample no. This is not the
576        preferred interface but is provided for use by C code.
577        The bufferg parameter is only required for LazyData.
578        \param sampleNo - Input - the given sample no.
579        \param bufferg - A buffer to compute (and store) sample data in will be selected from this group.
580        \return pointer to the sample data.
581    */
582      ESCRIPT_DLL_API
583      inline
584      const DataAbstract::ValueType::value_type*
585      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg=0);
586    
587    
588    /**    /**
589       \brief       \brief
590       Return the sample data for the given sample no. This is not the       Return the sample data for the given sample no. This is not the
591       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
592       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
593         \return pointer to the sample data.
594    */    */
595    ESCRIPT_DLL_API    ESCRIPT_DLL_API
596    inline    inline
597    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
598    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
599    {  
     return m_data->getSampleData(sampleNo);  
   }  
600    
601    /**    /**
602       \brief       \brief
# Line 521  class Data { Line 614  class Data {
614    
615    /**    /**
616       \brief       \brief
617       Return a view into the data for the data point specified.       Return a reference into the DataVector which points to the specified data point.
618       NOTE: Construction of the DataArrayView is a relatively expensive       \param sampleNo - Input -
619       operation.       \param dataPointNo - Input -
620      */
621      ESCRIPT_DLL_API
622      DataTypes::ValueType::const_reference
623      getDataPointRO(int sampleNo, int dataPointNo);
624    
625      /**
626         \brief
627         Return a reference into the DataVector which points to the specified data point.
628       \param sampleNo - Input -       \param sampleNo - Input -
629       \param dataPointNo - Input -       \param dataPointNo - Input -
630    */    */
631    ESCRIPT_DLL_API    ESCRIPT_DLL_API
632      DataTypes::ValueType::reference
633      getDataPointRW(int sampleNo, int dataPointNo);
634    
635    
636    
637      /**
638         \brief
639         Return the offset for the given sample and point within the sample
640      */
641      ESCRIPT_DLL_API
642    inline    inline
643    DataArrayView    DataTypes::ValueType::size_type
644    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
645                 int dataPointNo)                 int dataPointNo)
646    {    {
647                  return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
648    }    }
649    
650    /**    /**
# Line 541  class Data { Line 652  class Data {
652       Return a reference to the data point shape.       Return a reference to the data point shape.
653    */    */
654    ESCRIPT_DLL_API    ESCRIPT_DLL_API
655    const DataArrayView::ShapeType&    inline
656    getDataPointShape() const;    const DataTypes::ShapeType&
657      getDataPointShape() const
658      {
659        return m_data->getShape();
660      }
661    
662    /**    /**
663       \brief       \brief
# Line 566  class Data { Line 681  class Data {
681       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
682    */    */
683    ESCRIPT_DLL_API    ESCRIPT_DLL_API
684    DataArrayView::ValueType::size_type    DataTypes::ValueType::size_type
685    getLength() const;    getLength() const;
686    
687      /**
688      \brief Return true if this object contains no samples.
689      This is not the same as isEmpty()
690      */
691      ESCRIPT_DLL_API
692      bool
693      hasNoSamples() const
694      {
695        return getLength()==0;
696      }
697    
698    /**    /**
699       \brief       \brief
700       Assign the given value to the tag assocciated with name. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
701       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
702       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
703       \param tagKey - Input - Integer key.       \param name - Input - name of tag.
704       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
     ==>*  
705    */    */
706    ESCRIPT_DLL_API    ESCRIPT_DLL_API
707    void    void
# Line 605  class Data { Line 728  class Data {
728       object to type DataTagged if it is constant.       object to type DataTagged if it is constant.
729    
730       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
731         \param pointshape - Input - The shape of the value parameter
732       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
733      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
734    */    */
735    ESCRIPT_DLL_API    ESCRIPT_DLL_API
736    void    void
737    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
738                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
739                            const DataTypes::ValueType& value,
740                int dataOffset=0);
741    
742    
743    
744    /**    /**
745      \brief      \brief
# Line 644  class Data { Line 772  class Data {
772    ESCRIPT_DLL_API    ESCRIPT_DLL_API
773    Data    Data
774    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
775    
776    
777      ESCRIPT_DLL_API
778      Data
779      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
780                           double undef, Data& B, double Bmin, double Bstep,bool check_boundaries);
781    
782      ESCRIPT_DLL_API
783      Data
784      interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
785                           double undef,bool check_boundaries);
786    
787    
788    
789    
790      ESCRIPT_DLL_API
791      Data
792      interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
793                            Data& B, double Bmin, double Bstep, double undef,bool check_boundaries);
794    
795      ESCRIPT_DLL_API
796      Data
797      interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
798                            double undef,bool check_boundaries);
799    
800    /**    /**
801       \brief       \brief
802       Calculates the gradient of the data at the data points of functionspace.       Calculates the gradient of the data at the data points of functionspace.
# Line 659  class Data { Line 812  class Data {
812    grad() const;    grad() const;
813    
814    /**    /**
815       \brief      \brief
816       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
817       *    */
818      ESCRIPT_DLL_API
819      boost::python::object
820      integrateToTuple_const() const;
821    
822    
823      /**
824        \brief
825         Calculate the integral over the function space domain as a python tuple.
826    */    */
827    ESCRIPT_DLL_API    ESCRIPT_DLL_API
828    boost::python::numeric::array    boost::python::object
829    integrate() const;    integrateToTuple();
830    
831    
832    
833    /**    /**
834       \brief       \brief
# Line 732  class Data { Line 895  class Data {
895    /**    /**
896       \brief       \brief
897       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
898       *  
899         The method is not const because lazy data needs to be expanded before Lsup can be computed.
900         The _const form can be used when the Data object is const, however this will only work for
901         Data which is not Lazy.
902    
903         For Data which contain no samples (or tagged Data for which no tags in use have a value)
904         zero is returned.
905    */    */
906    ESCRIPT_DLL_API    ESCRIPT_DLL_API
907    double    double
908    Lsup() const;    Lsup();
909    
910      ESCRIPT_DLL_API
911      double
912      Lsup_const() const;
913    
914    
915    /**    /**
916       \brief       \brief
917       Return the maximum value of this Data object.       Return the maximum value of this Data object.
918       *  
919         The method is not const because lazy data needs to be expanded before sup can be computed.
920         The _const form can be used when the Data object is const, however this will only work for
921         Data which is not Lazy.
922    
923         For Data which contain no samples (or tagged Data for which no tags in use have a value)
924         a large negative value is returned.
925    */    */
926    ESCRIPT_DLL_API    ESCRIPT_DLL_API
927    double    double
928    sup() const;    sup();
929    
930      ESCRIPT_DLL_API
931      double
932      sup_const() const;
933    
934    
935    /**    /**
936       \brief       \brief
937       Return the minimum value of this Data object.       Return the minimum value of this Data object.
938       *  
939         The method is not const because lazy data needs to be expanded before inf can be computed.
940         The _const form can be used when the Data object is const, however this will only work for
941         Data which is not Lazy.
942    
943         For Data which contain no samples (or tagged Data for which no tags in use have a value)
944         a large positive value is returned.
945    */    */
946    ESCRIPT_DLL_API    ESCRIPT_DLL_API
947    double    double
948    inf() const;    inf();
949    
950      ESCRIPT_DLL_API
951      double
952      inf_const() const;
953    
954    
955    
956    /**    /**
957       \brief       \brief
# Line 786  class Data { Line 983  class Data {
983    /**    /**
984       \brief       \brief
985       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
986       the minimum value in this Data object.       the minimum component value in this Data object.
987         \note If you are working in python, please consider using Locator
988    instead of manually manipulating process and point IDs.
989    */    */
990    ESCRIPT_DLL_API    ESCRIPT_DLL_API
991    const boost::python::tuple    const boost::python::tuple
992    minGlobalDataPoint() const;    minGlobalDataPoint() const;
993    
994      /**
995         \brief
996         Return the (sample number, data-point number) of the data point with
997         the minimum component value in this Data object.
998         \note If you are working in python, please consider using Locator
999    instead of manually manipulating process and point IDs.
1000      */
1001    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1002    void    const boost::python::tuple
1003    calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;    maxGlobalDataPoint() const;
1004    
1005    
1006    
1007    /**    /**
1008       \brief       \brief
1009       Return the sign of each data point of this Data object.       Return the sign of each data point of this Data object.
# Line 1095  class Data { Line 1304  class Data {
1304    void    void
1305    saveVTK(std::string fileName) const;    saveVTK(std::string fileName) const;
1306    
1307    
1308    
1309    /**    /**
1310       \brief       \brief
1311       Overloaded operator +=       Overloaded operator +=
# Line 1212  class Data { Line 1423  class Data {
1423    */    */
1424    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1425    Data    Data
1426    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1427    
1428    /**    /**
1429       \brief       \brief
# Line 1225  class Data { Line 1436  class Data {
1436    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1437    void    void
1438    setSlice(const Data& value,    setSlice(const Data& value,
1439             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Archive the current Data object to the given file.  
      \param fileName - Input - file to archive to.  
   */  
   ESCRIPT_DLL_API  
   void  
   archiveData(const std::string fileName);  
   
   /**  
      \brief  
      Extract the Data object archived in the given file, overwriting  
      the current Data object.  
      Note - the current object must be of type DataEmpty.  
      \param fileName - Input - file to extract from.  
      \param fspace - Input - a suitable FunctionSpace descibing the data.  
   */  
   ESCRIPT_DLL_API  
   void  
   extractData(const std::string fileName,  
               const FunctionSpace& fspace);  
   
1440    
1441    /**    /**
1442       \brief       \brief
# Line 1290  class Data { Line 1478  class Data {
1478    /**    /**
1479       \brief       \brief
1480       return the object produced by the factory, which is a DataConstant or DataExpanded       return the object produced by the factory, which is a DataConstant or DataExpanded
1481        TODO Ownership of this object should be explained in doco.
1482    */    */
1483    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1484          DataAbstract*          DataAbstract*
1485          borrowData(void) const;          borrowData(void) const;
1486    
1487      ESCRIPT_DLL_API
1488            DataAbstract_ptr
1489            borrowDataPtr(void) const;
1490    
1491      ESCRIPT_DLL_API
1492            DataReady_ptr
1493            borrowReadyPtr(void) const;
1494    
1495    
1496    
1497      /**
1498         \brief
1499         Return a pointer to the beginning of the datapoint at the specified offset.
1500         TODO Eventually these should be inlined.
1501         \param i - position(offset) in the underlying datastructure
1502      */
1503    
1504      ESCRIPT_DLL_API
1505            DataTypes::ValueType::const_reference
1506            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1507    
1508    
1509      ESCRIPT_DLL_API
1510            DataTypes::ValueType::reference
1511            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1512    
1513    
1514    
1515    /**
1516       \brief Create a buffer for use by getSample
1517       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1518       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1519      
1520       In multi-threaded sections, this needs to be called on each thread.
1521    
1522       \return A BufferGroup* if Data is lazy, NULL otherwise.
1523       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1524    */
1525      ESCRIPT_DLL_API
1526      BufferGroup*
1527      allocSampleBuffer() const;
1528    
1529    /**
1530       \brief Free a buffer allocated with allocSampleBuffer.
1531       \param buffer Input - pointer to the buffer to deallocate.
1532    */
1533    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1534    
1535   protected:   protected:
1536    
1537   private:   private:
1538    
1539    template <class BinaryOp>
1540      double
1541      lazyAlgWorker(double init, int mpiop_type);
1542    
1543      double
1544      LsupWorker() const;
1545    
1546      double
1547      supWorker() const;
1548    
1549      double
1550      infWorker() const;
1551    
1552      boost::python::object
1553      integrateWorker() const;
1554    
1555      void
1556      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1557    
1558      void
1559      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1560    
1561    
1562    /**    /**
1563       \brief       \brief
1564       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1370  class Data { Line 1630  class Data {
1630       \brief       \brief
1631       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1632    */    */
1633    template <class IValueType>  
1634    void    void
1635    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1636             const DataTypes::ShapeType& shape,
1637               const FunctionSpace& what,               const FunctionSpace& what,
1638               bool expanded);               bool expanded);
1639    
1640      void
1641      initialise(const WrappedArray& value,
1642                     const FunctionSpace& what,
1643                     bool expanded);
1644    
1645    //    //
1646    // flag to protect the data object against any update    // flag to protect the data object against any update
1647    bool m_protected;    bool m_protected;
1648      mutable bool m_shared;
1649      bool m_lazy;
1650    
1651    //    //
1652    // pointer to the actual data object    // pointer to the actual data object
1653    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1654      DataAbstract_ptr m_data;
1655    
1656  };  // If possible please use getReadyPtr instead.
1657    // But see warning below.
1658      const DataReady*
1659      getReady() const;
1660    
1661  template <class IValueType>    DataReady*
1662  void    getReady();
1663  Data::initialise(const IValueType& value,  
1664                   const FunctionSpace& what,  
1665                   bool expanded)  // Be wary of using this for local operations since it (temporarily) increases reference count.
1666  {  // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1667    // getReady() instead
1668      DataReady_ptr
1669      getReadyPtr();
1670    
1671      const_DataReady_ptr
1672      getReadyPtr() const;
1673    
1674    
1675      /**
1676       \brief Update the Data's shared flag
1677       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1678       For internal use only.
1679      */
1680      void updateShareStatus(bool nowshared) const
1681      {
1682        m_shared=nowshared;     // m_shared is mutable
1683      }
1684    
1685      // In the isShared() method below:
1686      // A problem would occur if m_data (the address pointed to) were being modified
1687      // while the call m_data->is_shared is being executed.
1688      //
1689      // Q: So why do I think this code can be thread safe/correct?
1690      // A: We need to make some assumptions.
1691      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1692      // 2. We assume that no constructions or assignments which will share previously unshared
1693      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1694    //    //
1695    // Construct a Data object of the appropriate type.    // This means that the only transition we need to consider, is when a previously shared object is
1696    // Construct the object first as there seems to be a bug which causes    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1697    // undefined behaviour if an exception is thrown during construction    // In those cases the m_shared flag changes to false after m_data has completed changing.
1698    // within the shared_ptr constructor.    // For any threads executing before the flag switches they will assume the object is still shared.
1699    if (expanded) {    bool isShared() const
1700      DataAbstract* temp=new DataExpanded(value,what);    {
1701      boost::shared_ptr<DataAbstract> temp_data(temp);      return m_shared;
1702      m_data=temp_data;  /*  if (m_shared) return true;
1703    } else {      if (m_data->isShared())        
1704      DataAbstract* temp=new DataConstant(value,what);      {                  
1705      boost::shared_ptr<DataAbstract> temp_data(temp);          updateShareStatus(true);
1706      m_data=temp_data;          return true;
1707        }
1708        return false;*/
1709    }    }
1710    
1711      void forceResolve()
1712      {
1713        if (isLazy())
1714        {
1715            #ifdef _OPENMP
1716            if (omp_in_parallel())
1717            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1718            throw DataException("Please do not call forceResolve() in a parallel region.");
1719            }
1720            #endif
1721            resolve();
1722        }
1723      }
1724    
1725      /**
1726      \brief if another object is sharing out member data make a copy to work with instead.
1727      This code should only be called from single threaded sections of code.
1728      */
1729      void exclusiveWrite()
1730      {
1731    #ifdef _OPENMP
1732      if (omp_in_parallel())
1733      {
1734    // *((int*)0)=17;
1735        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1736      }
1737    #endif
1738        forceResolve();
1739        if (isShared())
1740        {
1741            DataAbstract* t=m_data->deepCopy();
1742            set_m_data(DataAbstract_ptr(t));
1743        }
1744      }
1745    
1746      /**
1747      \brief checks if caller can have exclusive write to the object
1748      */
1749      void checkExclusiveWrite()
1750      {
1751        if  (isLazy() || isShared())
1752        {
1753            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1754        }
1755      }
1756    
1757      /**
1758      \brief Modify the data abstract hosted by this Data object
1759      For internal use only.
1760      Passing a pointer to null is permitted (do this in the destructor)
1761      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1762      */
1763      void set_m_data(DataAbstract_ptr p);
1764    
1765      friend class DataAbstract;        // To allow calls to updateShareStatus
1766    
1767    };
1768    
1769    }   // end namespace escript
1770    
1771    
1772    // No, this is not supposed to be at the top of the file
1773    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1774    // so that I can dynamic cast between them below.
1775    #include "DataReady.h"
1776    #include "DataLazy.h"
1777    
1778    namespace escript
1779    {
1780    
1781    inline
1782    const DataReady*
1783    Data::getReady() const
1784    {
1785       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1786       EsysAssert((dr!=0), "Error - casting to DataReady.");
1787       return dr;
1788    }
1789    
1790    inline
1791    DataReady*
1792    Data::getReady()
1793    {
1794       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1795       EsysAssert((dr!=0), "Error - casting to DataReady.");
1796       return dr;
1797    }
1798    
1799    // Be wary of using this for local operations since it (temporarily) increases reference count.
1800    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1801    // getReady() instead
1802    inline
1803    DataReady_ptr
1804    Data::getReadyPtr()
1805    {
1806       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1807       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1808       return dr;
1809    }
1810    
1811    
1812    inline
1813    const_DataReady_ptr
1814    Data::getReadyPtr() const
1815    {
1816       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1817       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1818       return dr;
1819  }  }
1820    
1821    inline
1822    DataAbstract::ValueType::value_type*
1823    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1824    {
1825       if (isLazy())
1826       {
1827        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1828       }
1829       return getReady()->getSampleDataRW(sampleNo);
1830    }
1831    
1832    inline
1833    const DataAbstract::ValueType::value_type*
1834    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1835    {
1836       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1837       if (l!=0)
1838       {
1839        size_t offset=0;
1840        if (bufferg==NULL)
1841        {
1842            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1843        }
1844        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1845        return &((*res)[offset]);
1846       }
1847       return getReady()->getSampleDataRO(sampleNo);
1848    }
1849    
1850    
1851    
1852    /**
1853       Modify a filename for MPI parallel output to multiple files
1854    */
1855    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1856    
1857  /**  /**
1858     Binary Data object operators.     Binary Data object operators.
1859  */  */
# Line 1519  ESCRIPT_DLL_API std::ostream& operator<< Line 1965  ESCRIPT_DLL_API std::ostream& operator<<
1965  /**  /**
1966    \brief    \brief
1967    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1968    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1969    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1970    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1971    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1972  */  */
1973  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1974  Data  Data
1975  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1976                       Data& arg1,                       Data& arg_1,
1977                       int axis_offset=0,                       int axis_offset=0,
1978                       int transpose=0);                       int transpose=0);
1979    
   
 /**  
   \brief  
   Compute a tensor operation with two Data objects  
   \param arg0 - Input - Data object  
   \param arg1 - Input - Data object  
   \param operation - Input - Binary op functor  
 */  
 template <typename BinaryFunction>  
 ESCRIPT_DLL_API  
 Data  
 C_TensorBinaryOperation(Data const &arg0,  
                         Data const &arg1,  
                         BinaryFunction operation);  
   
 /**  
   \brief  
   Return true if operands are equivalent, else return false.  
   NB: this operator does very little at this point, and isn't to  
   be relied on. Requires further implementation.  
 */  
 // ESCRIPT_DLL_API bool operator==(const Data& left, const Data& right);  
   
1980  /**  /**
1981    \brief    \brief
1982    Perform the given binary operation with this and right as operands.    Perform the given binary operation with this and right as operands.
# Line 1567  Data::binaryOp(const Data& right, Line 1990  Data::binaryOp(const Data& right,
1990  {  {
1991     //     //
1992     // 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
1993     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1994       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1995     }     }
1996    
1997       if (isLazy() || right.isLazy())
1998       {
1999         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
2000       }
2001     //     //
2002     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
2003     Data tempRight(right);     Data tempRight(right);
2004    
2005     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
2006       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
2007         //         //
# Line 1582  Data::binaryOp(const Data& right, Line 2011  Data::binaryOp(const Data& right,
2011         //         //
2012         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2013         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2014         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2015           set_m_data(tempLeft.m_data);
2016       }       }
2017     }     }
2018     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1598  Data::binaryOp(const Data& right, Line 2028  Data::binaryOp(const Data& right,
2028       // of any data type       // of any data type
2029       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2030       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2031       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2032     } else if (isTagged()) {     } else if (isTagged()) {
2033       //       //
2034       // 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 1646  Data::algorithm(BinaryFunction operation Line 2076  Data::algorithm(BinaryFunction operation
2076      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2077      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2078      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2079      } else if (isEmpty()) {
2080        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2081      } else if (isLazy()) {
2082        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2083      } else {
2084        throw DataException("Error - Data encapsulates an unknown type.");
2085    }    }
   return 0;  
2086  }  }
2087    
2088  /**  /**
# Line 1663  inline Line 2098  inline
2098  Data  Data
2099  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2100  {  {
2101    if (isExpanded()) {    if (isEmpty()) {
2102      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2103      }
2104      else if (isExpanded()) {
2105        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2106      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2107      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2108      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2109      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2110      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2111      return result;      return result;
2112    } else if (isTagged()) {    }
2113      else if (isTagged()) {
2114      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());  
2115      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2116      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2117        defval[0]=0;
2118        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2119      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2120      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2121    } else if (isConstant()) {    }
2122      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2123        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2124      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2125      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2126      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2127      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2128      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2129      return result;      return result;
2130      } else if (isLazy()) {
2131        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2132      } else {
2133        throw DataException("Error - Data encapsulates an unknown type.");
2134    }    }
   Data falseRetVal; // to keep compiler quiet  
   return falseRetVal;  
2135  }  }
2136    
2137    /**
2138      \brief
2139      Compute a tensor operation with two Data objects
2140      \param arg_0 - Input - Data object
2141      \param arg_1 - Input - Data object
2142      \param operation - Input - Binary op functor
2143    */
2144  template <typename BinaryFunction>  template <typename BinaryFunction>
2145    inline
2146  Data  Data
2147  C_TensorBinaryOperation(Data const &arg_0,  C_TensorBinaryOperation(Data const &arg_0,
2148                          Data const &arg_1,                          Data const &arg_1,
2149                          BinaryFunction operation)                          BinaryFunction operation)
2150  {  {
2151      if (arg_0.isEmpty() || arg_1.isEmpty())
2152      {
2153         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2154      }
2155      if (arg_0.isLazy() || arg_1.isLazy())
2156      {
2157         throw DataException("Error - Operations not permitted on lazy data.");
2158      }
2159    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2160    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
2161    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
# Line 1730  C_TensorBinaryOperation(Data const &arg_ Line 2177  C_TensorBinaryOperation(Data const &arg_
2177    // Get rank and shape of inputs    // Get rank and shape of inputs
2178    int rank0 = arg_0_Z.getDataPointRank();    int rank0 = arg_0_Z.getDataPointRank();
2179    int rank1 = arg_1_Z.getDataPointRank();    int rank1 = arg_1_Z.getDataPointRank();
2180    DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();    DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2181    DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();    DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2182    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2183    int size1 = arg_1_Z.getDataPointSize();    int size1 = arg_1_Z.getDataPointSize();
   
2184    // Declare output Data object    // Declare output Data object
2185    Data res;    Data res;
2186    
2187    if (shape0 == shape1) {    if (shape0 == shape1) {
   
2188      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2189        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2190        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2191        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2192        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2193    
2194        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2195      }      }
2196      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
# Line 1762  C_TensorBinaryOperation(Data const &arg_ Line 2208  C_TensorBinaryOperation(Data const &arg_
2208    
2209        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2210        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2211        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2212        // Get the views  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2213        // Get the pointers to the actual data        // Get the pointers to the actual data
2214        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2215        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2216    
2217        // Compute a result for the default        // Compute a result for the default
2218        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2219        // Compute a result for each tag        // Compute a result for each tag
2220        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2221        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2222        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2223          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2224          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2225          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2226          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2227          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2228        }        }
2229    
2230      }      }
2231      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
   
2232        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2233        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2234        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
# Line 1795  C_TensorBinaryOperation(Data const &arg_ Line 2238  C_TensorBinaryOperation(Data const &arg_
2238        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2239        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2240        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2241          res.requireWrite();
2242        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2243        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2244          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2245            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2246            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2247            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2248            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2249            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2250            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2251          }          }
2252        }        }
2253    
2254      }      }
2255      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
   
2256        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2257        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2258    
# Line 1823  C_TensorBinaryOperation(Data const &arg_ Line 2266  C_TensorBinaryOperation(Data const &arg_
2266    
2267        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2268        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2269        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);  
2270        // Get the views        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2271        // Get the pointers to the actual data        // Get the pointers to the actual data
2272        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2273        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2274        // Compute a result for the default        // Compute a result for the default
2275        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2276        // Compute a result for each tag        // Compute a result for each tag
2277        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2278        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2279        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2280          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2281          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2282          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
         double *ptr_0 = &view_0.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2283          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2284        }        }
2285    
2286      }      }
2287      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
   
2288        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2289        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2290    
# Line 1858  C_TensorBinaryOperation(Data const &arg_ Line 2296  C_TensorBinaryOperation(Data const &arg_
2296        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2297        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2298    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2299        // Get the pointers to the actual data        // Get the pointers to the actual data
2300        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2301        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2302        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2303    
2304        // Compute a result for the default        // Compute a result for the default
2305        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2306        // Merge the tags        // Merge the tags
# Line 1873  C_TensorBinaryOperation(Data const &arg_ Line 2308  C_TensorBinaryOperation(Data const &arg_
2308        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2309        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2310        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2311          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape          tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2312        }        }
2313        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2314          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2315        }        }
2316        // Compute a result for each tag        // Compute a result for each tag
2317        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2318        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2319          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);  
2320          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2321          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2322          double *ptr_0 = &view_0.getData(0);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2323          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2324          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2325        }        }
2326    
2327      }      }
2328      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
   
2329        // After finding a common function space above the two inputs have the same numSamples and num DPPS        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2330        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2331        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
# Line 1902  C_TensorBinaryOperation(Data const &arg_ Line 2335  C_TensorBinaryOperation(Data const &arg_
2335        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2336        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2337        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2338          res.requireWrite();
2339        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2340        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2341          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2342          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2343          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2344            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2345            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2346            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2347            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2348            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2349          }          }
2350        }        }
2351    
2352      }      }
2353      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2354        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2355        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2356        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
# Line 1927  C_TensorBinaryOperation(Data const &arg_ Line 2360  C_TensorBinaryOperation(Data const &arg_
2360        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2361        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2362        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2363          res.requireWrite();
2364        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2365        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2366          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2367            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2368            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2369            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);  
2370            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2371            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2372              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2373    
2374    
2375            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2376          }          }
2377        }        }
2378    
2379      }      }
2380      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
   
2381        // After finding a common function space above the two inputs have the same numSamples and num DPPS        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2382        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2383        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
# Line 1951  C_TensorBinaryOperation(Data const &arg_ Line 2387  C_TensorBinaryOperation(Data const &arg_
2387        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2388        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2389        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2390          res.requireWrite();
2391        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2392        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2393          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2394          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2395          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2396            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2397            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2398            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2399            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2400            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2401          }          }
2402        }        }
2403    
2404      }      }
2405      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
   
2406        // After finding a common function space above the two inputs have the same numSamples and num DPPS        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2407        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2408        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
# Line 1976  C_TensorBinaryOperation(Data const &arg_ Line 2412  C_TensorBinaryOperation(Data const &arg_
2412        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2413        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2414        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2415          res.requireWrite();
2416        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2417        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2418          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        dataPointNo_0=0;
2419    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2420            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2421            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2422            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2423            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2424            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2425            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2426            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2427          }  //       }
2428        }        }
2429    
2430      }      }
# Line 1995  C_TensorBinaryOperation(Data const &arg_ Line 2433  C_TensorBinaryOperation(Data const &arg_
2433      }      }
2434    
2435    } else if (0 == rank0) {    } else if (0 == rank0) {
   
2436      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2437        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2438        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2439        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2440        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2441        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2442      }      }
2443      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
# Line 2018  C_TensorBinaryOperation(Data const &arg_ Line 2455  C_TensorBinaryOperation(Data const &arg_
2455    
2456        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2457        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2458        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2459        // Get the views  
2460        DataArrayView view_1 = tmp_1->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2461        DataArrayView view_2 = tmp_2->getDefaultValue();        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2462        // Get the pointers to the actual data  
       double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2463        // Compute a result for the default        // Compute a result for the default
2464        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2465        // Compute a result for each tag        // Compute a result for each tag
2466        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2467        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2468        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2469          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2470          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2471          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2472          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2473        }        }
2474    
# Line 2051  C_TensorBinaryOperation(Data const &arg_ Line 2484  C_TensorBinaryOperation(Data const &arg_
2484        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2485        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2486        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2487          res.requireWrite();
2488        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2489        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2490          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2491            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2492            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2493            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2494            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2495            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2496            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2497    
2498          }          }
# Line 2080  C_TensorBinaryOperation(Data const &arg_ Line 2514  C_TensorBinaryOperation(Data const &arg_
2514    
2515        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2516        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2517        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2518        // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2519        // Get the pointers to the actual data        // Get the pointers to the actual data
2520        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2521        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2522    
2523    
2524        // Compute a result for the default        // Compute a result for the default
2525        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2526        // Compute a result for each tag        // Compute a result for each tag
2527        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2528        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2529        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2530          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2531          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2532          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2533          double *ptr_0 = &view_0.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2534          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2535        }        }
2536    
# Line 2115  C_TensorBinaryOperation(Data const &arg_ Line 2548  C_TensorBinaryOperation(Data const &arg_
2548        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2549        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2550    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2551        // Get the pointers to the actual data        // Get the pointers to the actual data
2552        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2553        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2554        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2555    
2556        // Compute a result for the default        // Compute a result for the default
2557        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2558        // Merge the tags        // Merge the tags
# Line 2130  C_TensorBinaryOperation(Data const &arg_ Line 2560  C_TensorBinaryOperation(Data const &arg_
2560        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2561        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2562        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2563          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape          tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2564        }        }
2565        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2566          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2567        }        }
2568        // Compute a result for each tag        // Compute a result for each tag
2569        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2570        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2571          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2572          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2573          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2574          double *ptr_0 = &view_0.getData(0);  
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2575          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2576        }        }
2577    
# Line 2159  C_TensorBinaryOperation(Data const &arg_ Line 2587  C_TensorBinaryOperation(Data const &arg_
2587        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2588        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2589        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2590          res.requireWrite();
2591        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2592        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2593          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2594          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2595          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2596            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2597            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2598            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2599            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2600            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2601          }          }
2602        }        }
2603    
2604      }      }
2605      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2606        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2607        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2608        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
# Line 2184  C_TensorBinaryOperation(Data const &arg_ Line 2612  C_TensorBinaryOperation(Data const &arg_
2612        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2613        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2614        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2615          res.requireWrite();
2616        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2617        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2618          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2619            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2620            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2621            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2622            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2623            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2624            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2625          }          }
2626        }        }
# Line 2209  C_TensorBinaryOperation(Data const &arg_ Line 2638  C_TensorBinaryOperation(Data const &arg_
2638        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2639        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2640        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2641          res.requireWrite();
2642        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2643        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2644          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2645          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2646          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2647            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2648            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2649            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2650            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2651            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2652          }          }
2653        }        }
# Line 2234  C_TensorBinaryOperation(Data const &arg_ Line 2664  C_TensorBinaryOperation(Data const &arg_
2664        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2665        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2666        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2667          res.requireWrite();
2668        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2669        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2670          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2671            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2672            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2673            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2674            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2675            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2676            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2677            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2678          }          }
2679        }        }
# Line 2253  C_TensorBinaryOperation(Data const &arg_ Line 2684  C_TensorBinaryOperation(Data const &arg_
2684      }      }
2685    
2686    } else if (0 == rank1) {    } else if (0 == rank1) {
   
2687      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2688        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2689        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2690        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2691        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2692        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2693      }      }
2694      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
# Line 2276  C_TensorBinaryOperation(Data const &arg_ Line 2706  C_TensorBinaryOperation(Data const &arg_
2706    
2707        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2708        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2709        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2710        // Get the views  
2711        DataArrayView view_1 = tmp_1->getDefaultValue();        //Get the pointers to the actual data
2712        DataArrayView view_2 = tmp_2->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2713        // Get the pointers to the actual data        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2714        double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2715        // Compute a result for the default        // Compute a result for the default
2716        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2717        // Compute a result for each tag        // Compute a result for each tag
2718        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2719        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2720        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2721          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2722          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2723          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2724          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2725        }        }
   
2726      }      }
2727      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2728    
# Line 2309  C_TensorBinaryOperation(Data const &arg_ Line 2735  C_TensorBinaryOperation(Data const &arg_
2735        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2736        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2737        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2738          res.requireWrite();
2739        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2740        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2741          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2742            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2743            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2744            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2745            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2746            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2747            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2748          }          }
2749        }        }
# Line 2337  C_TensorBinaryOperation(Data const &arg_ Line 2764  C_TensorBinaryOperation(Data const &arg_
2764    
2765        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2766        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2767        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2768        // Get the pointers to the actual data        // Get the pointers to the actual data
2769        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2770        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2771        // Compute a result for the default        // Compute a result for the default
2772        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2773        // Compute a result for each tag        // Compute a result for each tag
2774        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2775        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2776        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2777          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2778          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2779          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
         double *ptr_0 = &view_0.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2780          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2781        }        }
2782    
# Line 2372  C_TensorBinaryOperation(Data const &arg_ Line 2794  C_TensorBinaryOperation(Data const &arg_
2794        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2795        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2796    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2797        // Get the pointers to the actual data        // Get the pointers to the actual data
2798        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2799        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2800        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2801    
2802        // Compute a result for the default        // Compute a result for the default
2803        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2804        // Merge the tags        // Merge the tags
# Line 2387  C_TensorBinaryOperation(Data const &arg_ Line 2806  C_TensorBinaryOperation(Data const &arg_
2806        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2807        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2808        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2809          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape          tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2810        }        }
2811        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2812          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2813        }        }
2814        // Compute a result for each tag        // Compute a result for each tag
2815        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2816        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2817          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2818          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2819          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
         double *ptr_0 = &view_0.getData(0);  
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2820          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2821        }        }
2822    
# Line 2416  C_TensorBinaryOperation(Data const &arg_ Line 2832  C_TensorBinaryOperation(Data const &arg_
2832        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2833        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2834        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2835          res.requireWrite();
2836        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2837        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2838          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2839          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2840          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2841            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2842            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2843            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2844            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2845            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2846          }          }
2847        }        }
2848    
2849      }      }
2850      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2851        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2852        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2853        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
# Line 2441  C_TensorBinaryOperation(Data const &arg_ Line 2857  C_TensorBinaryOperation(Data const &arg_
2857        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2858        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2859        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2860          res.requireWrite();
2861        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2862        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2863          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2864            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2865            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2866            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2867            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2868            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2869            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2870          }          }
2871        }        }
# Line 2466  C_TensorBinaryOperation(Data const &arg_ Line 2883  C_TensorBinaryOperation(Data const &arg_
2883        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2884        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2885        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2886          res.requireWrite();
2887        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2888        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2889          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2890          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2891          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2892            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2893            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2894            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2895            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2896            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2897          }          }
2898        }        }
# Line 2491  C_TensorBinaryOperation(Data const &arg_ Line 2909  C_TensorBinaryOperation(Data const &arg_
2909        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2910        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2911        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2912          res.requireWrite();
2913        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2914        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2915          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2916            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2917            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2918            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2919            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2920            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2921            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2922            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2923          }          }
2924        }        }
# Line 2521  Data Line 2940  Data
2940  C_TensorUnaryOperation(Data const &arg_0,  C_TensorUnaryOperation(Data const &arg_0,
2941                         UnaryFunction operation)                         UnaryFunction operation)
2942  {  {
2943      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2944      {
2945         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2946      }
2947      if (arg_0.isLazy())
2948      {
2949         throw DataException("Error - Operations not permitted on lazy data.");
2950      }
2951    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2952    Data arg_0_Z = Data(arg_0);    Data arg_0_Z = Data(arg_0);
2953    
2954    // Get rank and shape of inputs    // Get rank and shape of inputs
2955    int rank0 = arg_0_Z.getDataPointRank();    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
   DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();  
2956    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2957    
2958    // Declare output Data object    // Declare output Data object
# Line 2534  C_TensorUnaryOperation(Data const &arg_0 Line 2960  C_TensorUnaryOperation(Data const &arg_0
2960    
2961    if (arg_0_Z.isConstant()) {    if (arg_0_Z.isConstant()) {
2962      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2963      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2964      double *ptr_2 = &((res.getPointDataView().getData())[0]);      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2965      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2966    }    }
2967    else if (arg_0_Z.isTagged()) {    else if (arg_0_Z.isTagged()) {
# Line 2548  C_TensorUnaryOperation(Data const &arg_0 Line 2974  C_TensorUnaryOperation(Data const &arg_0
2974      res.tag();      res.tag();
2975      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2976    
     // Get the views  
     DataArrayView view_0 = tmp_0->getDefaultValue();  
     DataArrayView view_2 = tmp_2->getDefaultValue();  
2977      // Get the pointers to the actual data      // Get the pointers to the actual data
2978      double *ptr_0 = &((view_0.getData())[0]);      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2979      double *ptr_2 = &((view_2.getData())[0]);      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2980      // Compute a result for the default      // Compute a result for the default
2981      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2982      // Compute a result for each tag      // Compute a result for each tag
2983      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2984      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2985      for (i=lookup_0.begin();i!=lookup_0.end();i++) {      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2986        tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());        tmp_2->addTag(i->first);
2987        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2988        DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
       double *ptr_0 = &view_0.getData(0);  
       double *ptr_2 = &view_2.getData(0);  
2989        tensor_unary_operation(size0, ptr_0, ptr_2, operation);        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2990      }      }
2991    
# Line 2580  C_TensorUnaryOperation(Data const &arg_0 Line 3001  C_TensorUnaryOperation(Data const &arg_0
3001      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3002      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3003      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3004        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {      dataPointNo_0=0;
3005    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3006          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3007          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3008          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3009          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3010          tensor_unary_operation(size0, ptr_0, ptr_2, operation);          tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3011        }  //      }
3012      }      }
   
3013    }    }
3014    else {    else {
3015      throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");      throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");

Legend:
Removed from v.1387  
changed lines
  Added in v.2721

  ViewVC Help
Powered by ViewVC 1.1.26