/[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 2742 by jfenwick, Thu Nov 12 06:03:37 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);
   
   /**  
      \brief  
      Constructor which will create Tagged data if expanded is false.  
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
     ==>*  
   */  
   ESCRIPT_DLL_API  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   ESCRIPT_DLL_API  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
145    
146    /**    /**
147       \brief       \brief
148       Constructor which copies data from any object that can be converted into       Constructor which copies data from any object that can be treated like a python array/sequence.
      a python numarray.  
149    
150       \param value - Input - Input data.       \param value - Input - Input data.
151       \param what - Input - A description of what this data represents.       \param what - Input - A description of what this data represents.
# Line 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.
      NOTE: Construction of the DataArrayView is a relatively expensive  
      operation.  
618       \param sampleNo - Input -       \param sampleNo - Input -
619       \param dataPointNo - Input -       \param dataPointNo - Input -
620    */    */
621    ESCRIPT_DLL_API    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 -
629         \param dataPointNo - Input -
630      */
631      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    ESCRIPT_DLL_API
819    boost::python::numeric::array    boost::python::object
820    integrate() const;    integrateToTuple_const() const;
821    
822    
823      /**
824        \brief
825         Calculate the integral over the function space domain as a python tuple.
826      */
827      ESCRIPT_DLL_API
828      boost::python::object
829      integrateToTuple();
830    
831    
832    
833    /**    /**
834       \brief       \brief
# Line 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 1143  class Data { Line 1354  class Data {
1354    Data& operator/=(const boost::python::object& right);    Data& operator/=(const boost::python::object& right);
1355    
1356    /**    /**
1357        \brief return inverse of matricies.
1358      */
1359      ESCRIPT_DLL_API
1360      Data
1361      matrixInverse() const;
1362    
1363      /**
1364       \brief       \brief
1365       Returns true if this can be interpolated to functionspace.       Returns true if this can be interpolated to functionspace.
1366    */    */
# Line 1212  class Data { Line 1430  class Data {
1430    */    */
1431    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1432    Data    Data
1433    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1434    
1435    /**    /**
1436       \brief       \brief
# Line 1225  class Data { Line 1443  class Data {
1443    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1444    void    void
1445    setSlice(const Data& value,    setSlice(const Data& value,
1446             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
   
   /**  
      \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);  
   
1447    
1448    /**    /**
1449       \brief       \brief
# Line 1290  class Data { Line 1485  class Data {
1485    /**    /**
1486       \brief       \brief
1487       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
1488        TODO Ownership of this object should be explained in doco.
1489    */    */
1490    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1491          DataAbstract*          DataAbstract*
1492          borrowData(void) const;          borrowData(void) const;
1493    
1494      ESCRIPT_DLL_API
1495            DataAbstract_ptr
1496            borrowDataPtr(void) const;
1497    
1498      ESCRIPT_DLL_API
1499            DataReady_ptr
1500            borrowReadyPtr(void) const;
1501    
1502    
1503    
1504      /**
1505         \brief
1506         Return a pointer to the beginning of the datapoint at the specified offset.
1507         TODO Eventually these should be inlined.
1508         \param i - position(offset) in the underlying datastructure
1509      */
1510    
1511      ESCRIPT_DLL_API
1512            DataTypes::ValueType::const_reference
1513            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1514    
1515    
1516      ESCRIPT_DLL_API
1517            DataTypes::ValueType::reference
1518            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1519    
1520    
1521    
1522    /**
1523       \brief Create a buffer for use by getSample
1524       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1525       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1526      
1527       In multi-threaded sections, this needs to be called on each thread.
1528    
1529       \return A BufferGroup* if Data is lazy, NULL otherwise.
1530       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1531    */
1532      ESCRIPT_DLL_API
1533      BufferGroup*
1534      allocSampleBuffer() const;
1535    
1536    /**
1537       \brief Free a buffer allocated with allocSampleBuffer.
1538       \param buffer Input - pointer to the buffer to deallocate.
1539    */
1540    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1541    
1542   protected:   protected:
1543    
1544   private:   private:
1545    
1546    template <class BinaryOp>
1547      double
1548    #ifdef PASO_MPI
1549      lazyAlgWorker(double init, MPI_Op mpiop_type);
1550    #else
1551      lazyAlgWorker(double init);
1552    #endif
1553    
1554      double
1555      LsupWorker() const;
1556    
1557      double
1558      supWorker() const;
1559    
1560      double
1561      infWorker() const;
1562    
1563      boost::python::object
1564      integrateWorker() const;
1565    
1566      void
1567      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1568    
1569      void
1570      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1571    
1572      // For internal use in Data.cpp only!
1573      // other uses should call the main entry points and allow laziness
1574      Data
1575      minval_nonlazy() const;
1576    
1577      // For internal use in Data.cpp only!
1578      Data
1579      maxval_nonlazy() const;
1580    
1581    
1582    /**    /**
1583       \brief       \brief
1584       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1370  class Data { Line 1650  class Data {
1650       \brief       \brief
1651       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1652    */    */
1653    template <class IValueType>  
1654    void    void
1655    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1656             const DataTypes::ShapeType& shape,
1657               const FunctionSpace& what,               const FunctionSpace& what,
1658               bool expanded);               bool expanded);
1659    
1660      void
1661      initialise(const WrappedArray& value,
1662                     const FunctionSpace& what,
1663                     bool expanded);
1664    
1665    //    //
1666    // flag to protect the data object against any update    // flag to protect the data object against any update
1667    bool m_protected;    bool m_protected;
1668      mutable bool m_shared;
1669      bool m_lazy;
1670    
1671    //    //
1672    // pointer to the actual data object    // pointer to the actual data object
1673    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1674      DataAbstract_ptr m_data;
1675    
1676  };  // If possible please use getReadyPtr instead.
1677    // But see warning below.
1678      const DataReady*
1679      getReady() const;
1680    
1681  template <class IValueType>    DataReady*
1682  void    getReady();
1683  Data::initialise(const IValueType& value,  
1684                   const FunctionSpace& what,  
1685                   bool expanded)  // Be wary of using this for local operations since it (temporarily) increases reference count.
1686  {  // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1687    // getReady() instead
1688      DataReady_ptr
1689      getReadyPtr();
1690    
1691      const_DataReady_ptr
1692      getReadyPtr() const;
1693    
1694    
1695      /**
1696       \brief Update the Data's shared flag
1697       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1698       For internal use only.
1699      */
1700      void updateShareStatus(bool nowshared) const
1701      {
1702        m_shared=nowshared;     // m_shared is mutable
1703      }
1704    
1705      // In the isShared() method below:
1706      // A problem would occur if m_data (the address pointed to) were being modified
1707      // while the call m_data->is_shared is being executed.
1708      //
1709      // Q: So why do I think this code can be thread safe/correct?
1710      // A: We need to make some assumptions.
1711      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1712      // 2. We assume that no constructions or assignments which will share previously unshared
1713      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1714    //    //
1715    // 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
1716    // 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.
1717    // 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.
1718    // within the shared_ptr constructor.    // For any threads executing before the flag switches they will assume the object is still shared.
1719    if (expanded) {    bool isShared() const
1720      DataAbstract* temp=new DataExpanded(value,what);    {
1721      boost::shared_ptr<DataAbstract> temp_data(temp);      return m_shared;
1722      m_data=temp_data;  /*  if (m_shared) return true;
1723    } else {      if (m_data->isShared())        
1724      DataAbstract* temp=new DataConstant(value,what);      {                  
1725      boost::shared_ptr<DataAbstract> temp_data(temp);          updateShareStatus(true);
1726      m_data=temp_data;          return true;
1727        }
1728        return false;*/
1729      }
1730    
1731      void forceResolve()
1732      {
1733        if (isLazy())
1734        {
1735            #ifdef _OPENMP
1736            if (omp_in_parallel())
1737            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1738            throw DataException("Please do not call forceResolve() in a parallel region.");
1739            }
1740            #endif
1741            resolve();
1742        }
1743      }
1744    
1745      /**
1746      \brief if another object is sharing out member data make a copy to work with instead.
1747      This code should only be called from single threaded sections of code.
1748      */
1749      void exclusiveWrite()
1750      {
1751    #ifdef _OPENMP
1752      if (omp_in_parallel())
1753      {
1754    // *((int*)0)=17;
1755        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1756      }
1757    #endif
1758        forceResolve();
1759        if (isShared())
1760        {
1761            DataAbstract* t=m_data->deepCopy();
1762            set_m_data(DataAbstract_ptr(t));
1763        }
1764      }
1765    
1766      /**
1767      \brief checks if caller can have exclusive write to the object
1768      */
1769      void checkExclusiveWrite()
1770      {
1771        if  (isLazy() || isShared())
1772        {
1773            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1774        }
1775    }    }
1776    
1777      /**
1778      \brief Modify the data abstract hosted by this Data object
1779      For internal use only.
1780      Passing a pointer to null is permitted (do this in the destructor)
1781      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1782      */
1783      void set_m_data(DataAbstract_ptr p);
1784    
1785      friend class DataAbstract;        // To allow calls to updateShareStatus
1786    
1787    };
1788    
1789    }   // end namespace escript
1790    
1791    
1792    // No, this is not supposed to be at the top of the file
1793    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1794    // so that I can dynamic cast between them below.
1795    #include "DataReady.h"
1796    #include "DataLazy.h"
1797    
1798    namespace escript
1799    {
1800    
1801    inline
1802    const DataReady*
1803    Data::getReady() const
1804    {
1805       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1806       EsysAssert((dr!=0), "Error - casting to DataReady.");
1807       return dr;
1808    }
1809    
1810    inline
1811    DataReady*
1812    Data::getReady()
1813    {
1814       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1815       EsysAssert((dr!=0), "Error - casting to DataReady.");
1816       return dr;
1817  }  }
1818    
1819    // Be wary of using this for local operations since it (temporarily) increases reference count.
1820    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1821    // getReady() instead
1822    inline
1823    DataReady_ptr
1824    Data::getReadyPtr()
1825    {
1826       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1827       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1828       return dr;
1829    }
1830    
1831    
1832    inline
1833    const_DataReady_ptr
1834    Data::getReadyPtr() const
1835    {
1836       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1837       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1838       return dr;
1839    }
1840    
1841    inline
1842    DataAbstract::ValueType::value_type*
1843    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1844    {
1845       if (isLazy())
1846       {
1847        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1848       }
1849       return getReady()->getSampleDataRW(sampleNo);
1850    }
1851    
1852    inline
1853    const DataAbstract::ValueType::value_type*
1854    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1855    {
1856       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1857       if (l!=0)
1858       {
1859        size_t offset=0;
1860        if (bufferg==NULL)
1861        {
1862            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1863        }
1864        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1865        return &((*res)[offset]);
1866       }
1867       return getReady()->getSampleDataRO(sampleNo);
1868    }
1869    
1870    
1871    
1872    /**
1873       Modify a filename for MPI parallel output to multiple files
1874    */
1875    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1876    
1877  /**  /**
1878     Binary Data object operators.     Binary Data object operators.
1879  */  */
# Line 1519  ESCRIPT_DLL_API std::ostream& operator<< Line 1985  ESCRIPT_DLL_API std::ostream& operator<<
1985  /**  /**
1986    \brief    \brief
1987    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1988    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1989    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1990    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1991    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1992  */  */
1993  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1994  Data  Data
1995  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1996                       Data& arg1,                       Data& arg_1,
1997                       int axis_offset=0,                       int axis_offset=0,
1998                       int transpose=0);                       int transpose=0);
1999    
   
 /**  
   \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);  
   
2000  /**  /**
2001    \brief    \brief
2002    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 2010  Data::binaryOp(const Data& right,
2010  {  {
2011     //     //
2012     // if this has a rank of zero promote it to the rank of the RHS     // if this has a rank of zero promote it to the rank of the RHS
2013     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
2014       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
2015     }     }
2016    
2017       if (isLazy() || right.isLazy())
2018       {
2019         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
2020       }
2021     //     //
2022     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
2023     Data tempRight(right);     Data tempRight(right);
2024    
2025     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
2026       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
2027         //         //
# Line 1582  Data::binaryOp(const Data& right, Line 2031  Data::binaryOp(const Data& right,
2031         //         //
2032         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2033         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2034         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2035           set_m_data(tempLeft.m_data);
2036       }       }
2037     }     }
2038     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1598  Data::binaryOp(const Data& right, Line 2048  Data::binaryOp(const Data& right,
2048       // of any data type       // of any data type
2049       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2050       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2051       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2052     } else if (isTagged()) {     } else if (isTagged()) {
2053       //       //
2054       // Tagged data is operated on serially, the right hand side can be       // Tagged data is operated on serially, the right hand side can be
# Line 1646  Data::algorithm(BinaryFunction operation Line 2096  Data::algorithm(BinaryFunction operation
2096      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2097      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2098      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2099      } else if (isEmpty()) {
2100        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2101      } else if (isLazy()) {
2102        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2103      } else {
2104        throw DataException("Error - Data encapsulates an unknown type.");
2105    }    }
   return 0;  
2106  }  }
2107    
2108  /**  /**
# Line 1663  inline Line 2118  inline
2118  Data  Data
2119  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2120  {  {
2121    if (isExpanded()) {    if (isEmpty()) {
2122      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2123      }
2124      else if (isExpanded()) {
2125        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2126      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2127      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2128      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2129      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2130      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2131      return result;      return result;
2132    } else if (isTagged()) {    }
2133      else if (isTagged()) {
2134      DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
     DataArrayView::ShapeType viewShape;  
     DataArrayView::ValueType viewData(1);  
     viewData[0]=0;  
     DataArrayView defaultValue(viewData,viewShape);  
     DataTagged::TagListType keys;  
     DataTagged::ValueListType values;  
     DataTagged::DataMapType::const_iterator i;  
     for (i=dataT->getTagLookup().begin();i!=dataT->getTagLookup().end();i++) {  
       keys.push_back(i->first);  
       values.push_back(defaultValue);  
     }  
     Data result(keys,values,defaultValue,getFunctionSpace());  
     DataTagged* resultT=dynamic_cast<DataTagged*>(result.m_data.get());  
2135      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2136      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2137        defval[0]=0;
2138        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2139      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2140      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2141    } else if (isConstant()) {    }
2142      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2143        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2144      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2145      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2146      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2147      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2148      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2149      return result;      return result;
2150      } else if (isLazy()) {
2151        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2152      } else {
2153        throw DataException("Error - Data encapsulates an unknown type.");
2154    }    }
   Data falseRetVal; // to keep compiler quiet  
   return falseRetVal;  
2155  }  }
2156    
2157    /**
2158      \brief
2159      Compute a tensor operation with two Data objects
2160      \param arg_0 - Input - Data object
2161      \param arg_1 - Input - Data object
2162      \param operation - Input - Binary op functor
2163    */
2164  template <typename BinaryFunction>  template <typename BinaryFunction>
2165    inline
2166  Data  Data
2167  C_TensorBinaryOperation(Data const &arg_0,  C_TensorBinaryOperation(Data const &arg_0,
2168                          Data const &arg_1,                          Data const &arg_1,
2169                          BinaryFunction operation)                          BinaryFunction operation)
2170  {  {
2171      if (arg_0.isEmpty() || arg_1.isEmpty())
2172      {
2173         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2174      }
2175      if (arg_0.isLazy() || arg_1.isLazy())
2176      {
2177         throw DataException("Error - Operations not permitted on lazy data.");
2178      }
2179    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2180    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
2181    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
# Line 1730  C_TensorBinaryOperation(Data const &arg_ Line 2197  C_TensorBinaryOperation(Data const &arg_
2197    // Get rank and shape of inputs    // Get rank and shape of inputs
2198    int rank0 = arg_0_Z.getDataPointRank();    int rank0 = arg_0_Z.getDataPointRank();
2199    int rank1 = arg_1_Z.getDataPointRank();    int rank1 = arg_1_Z.getDataPointRank();
2200    DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();    DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2201    DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();    DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2202    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2203    int size1 = arg_1_Z.getDataPointSize();    int size1 = arg_1_Z.getDataPointSize();
   
2204    // Declare output Data object    // Declare output Data object
2205    Data res;    Data res;
2206    
2207    if (shape0 == shape1) {    if (shape0 == shape1) {
   
2208      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2209        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2210        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2211        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2212        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2213    
2214        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2215      }      }
2216      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 2228  C_TensorBinaryOperation(Data const &arg_
2228    
2229        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2230        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2231        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2232        // Get the views  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2233        // Get the pointers to the actual data        // Get the pointers to the actual data
2234        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2235        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2236    
2237        // Compute a result for the default        // Compute a result for the default
2238        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2239        // Compute a result for each tag        // Compute a result for each tag
2240        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2241        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2242        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2243          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2244          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2245          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2246          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2247          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2248        }        }
2249    
2250      }      }
2251      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
   
2252        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2253        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2254        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 2258  C_TensorBinaryOperation(Data const &arg_
2258        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2259        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2260        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2261          res.requireWrite();
2262        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2263        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2264          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2265            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2266            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2267            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2268            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2269            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2270            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2271          }          }
2272        }        }
2273    
2274      }      }
2275      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
   
2276        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2277        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2278    
# Line 1823  C_TensorBinaryOperation(Data const &arg_ Line 2286  C_TensorBinaryOperation(Data const &arg_
2286    
2287        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2288        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2289        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);  
2290        // 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();  
2291        // Get the pointers to the actual data        // Get the pointers to the actual data
2292        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2293        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2294        // Compute a result for the default        // Compute a result for the default
2295        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2296        // Compute a result for each tag        // Compute a result for each tag
2297        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2298        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2299        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2300          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2301          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2302          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);  
2303          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2304        }        }
2305    
2306      }      }
2307      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
   
2308        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2309        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2310    
# Line 1858  C_TensorBinaryOperation(Data const &arg_ Line 2316  C_TensorBinaryOperation(Data const &arg_
2316        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2317        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2318    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2319        // Get the pointers to the actual data        // Get the pointers to the actual data
2320        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2321        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2322        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2323    
2324        // Compute a result for the default        // Compute a result for the default
2325        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2326        // Merge the tags        // Merge the tags
# Line 1873  C_TensorBinaryOperation(Data const &arg_ Line 2328  C_TensorBinaryOperation(Data const &arg_
2328        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2329        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2330        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2331          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
2332        }        }
2333        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2334          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2335        }        }
2336        // Compute a result for each tag        // Compute a result for each tag
2337        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2338        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2339          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);  
2340          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2341          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2342          double *ptr_0 = &view_0.getData(0);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2343          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2344          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2345        }        }
2346    
2347      }      }
2348      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
   
2349        // After finding a common function space above the two inputs have the same numSamples and num DPPS        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2350        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2351        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 2355  C_TensorBinaryOperation(Data const &arg_
2355        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2356        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2357        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2358          res.requireWrite();
2359        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2360        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2361          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2362          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2363          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2364            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2365            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2366            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2367            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2368            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2369          }          }
2370        }        }
2371    
2372      }      }
2373      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2374        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2375        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2376        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 2380  C_TensorBinaryOperation(Data const &arg_
2380        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2381        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2382        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2383          res.requireWrite();
2384        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2385        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2386          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2387            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2388            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2389            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);  
2390            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2391            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2392              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2393    
2394    
2395            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2396          }          }
2397        }        }
2398    
2399      }      }
2400      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
   
2401        // After finding a common function space above the two inputs have the same numSamples and num DPPS        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2402        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2403        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 2407  C_TensorBinaryOperation(Data const &arg_
2407        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2408        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2409        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2410          res.requireWrite();
2411        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2412        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2413          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2414          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2415          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2416            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2417            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2418            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2419            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2420            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2421          }          }
2422        }        }
2423    
2424      }      }
2425      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
   
2426        // After finding a common function space above the two inputs have the same numSamples and num DPPS        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2427        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2428        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 2432  C_TensorBinaryOperation(Data const &arg_
2432        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2433        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2434        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2435          res.requireWrite();
2436        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2437        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2438          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        dataPointNo_0=0;
2439    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2440            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2441            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2442            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2443            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2444            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2445            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2446            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2447          }  //       }
2448        }        }
2449    
2450      }      }
# Line 1995  C_TensorBinaryOperation(Data const &arg_ Line 2453  C_TensorBinaryOperation(Data const &arg_
2453      }      }
2454    
2455    } else if (0 == rank0) {    } else if (0 == rank0) {
   
2456      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2457        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2458        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2459        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2460        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2461        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2462      }      }
2463      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 2475  C_TensorBinaryOperation(Data const &arg_
2475    
2476        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2477        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2478        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2479        // Get the views  
2480        DataArrayView view_1 = tmp_1->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2481        DataArrayView view_2 = tmp_2->getDefaultValue();        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2482        // Get the pointers to the actual data  
       double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2483        // Compute a result for the default        // Compute a result for the default
2484        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2485        // Compute a result for each tag        // Compute a result for each tag
2486        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2487        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2488        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2489          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2490          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2491          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);  
2492          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2493        }        }
2494    
# Line 2051  C_TensorBinaryOperation(Data const &arg_ Line 2504  C_TensorBinaryOperation(Data const &arg_
2504        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2505        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2506        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2507          res.requireWrite();
2508        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2509        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2510          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2511            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2512            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2513            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2514            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2515            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2516            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2517    
2518          }          }
# Line 2080  C_TensorBinaryOperation(Data const &arg_ Line 2534  C_TensorBinaryOperation(Data const &arg_
2534    
2535        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2536        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2537        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2538        // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2539        // Get the pointers to the actual data        // Get the pointers to the actual data
2540        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2541        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2542    
2543    
2544        // Compute a result for the default        // Compute a result for the default
2545        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2546        // Compute a result for each tag        // Compute a result for each tag
2547        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2548        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2549        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2550          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2551          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2552          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2553          double *ptr_0 = &view_0.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2554          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2555        }        }
2556    
# Line 2115  C_TensorBinaryOperation(Data const &arg_ Line 2568  C_TensorBinaryOperation(Data const &arg_
2568        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2569        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2570    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2571        // Get the pointers to the actual data        // Get the pointers to the actual data
2572        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2573        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2574        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2575    
2576        // Compute a result for the default        // Compute a result for the default
2577        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2578        // Merge the tags        // Merge the tags
# Line 2130  C_TensorBinaryOperation(Data const &arg_ Line 2580  C_TensorBinaryOperation(Data const &arg_
2580        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2581        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2582        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2583          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
2584        }        }
2585        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2586          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2587        }        }
2588        // Compute a result for each tag        // Compute a result for each tag
2589        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2590        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2591          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2592          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2593          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2594          double *ptr_0 = &view_0.getData(0);  
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2595          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2596        }        }
2597    
# Line 2159  C_TensorBinaryOperation(Data const &arg_ Line 2607  C_TensorBinaryOperation(Data const &arg_
2607        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2608        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2609        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2610          res.requireWrite();
2611        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2612        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2613          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2614          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2615          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2616            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2617            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2618            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2619            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2620            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2621          }          }
2622        }        }
2623    
2624      }      }
2625      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2626        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2627        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2628        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 2632  C_TensorBinaryOperation(Data const &arg_
2632        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2633        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2634        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2635          res.requireWrite();
2636        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2637        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2638          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2639            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2640            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2641            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2642            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2643            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2644            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2645          }          }
2646        }        }
# Line 2209  C_TensorBinaryOperation(Data const &arg_ Line 2658  C_TensorBinaryOperation(Data const &arg_
2658        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2659        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2660        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2661          res.requireWrite();
2662        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2663        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2664          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2665          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2666          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2667            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2668            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2669            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2670            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2671            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2672          }          }
2673        }        }
# Line 2234  C_TensorBinaryOperation(Data const &arg_ Line 2684  C_TensorBinaryOperation(Data const &arg_
2684        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2685        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2686        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2687          res.requireWrite();
2688        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2689        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2690          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2691            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2692            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2693            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2694            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2695            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2696            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2697            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2698          }          }
2699        }        }
# Line 2253  C_TensorBinaryOperation(Data const &arg_ Line 2704  C_TensorBinaryOperation(Data const &arg_
2704      }      }
2705    
2706    } else if (0 == rank1) {    } else if (0 == rank1) {
   
2707      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2708        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2709        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2710        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2711        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2712        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2713      }      }
2714      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 2726  C_TensorBinaryOperation(Data const &arg_
2726    
2727        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2728        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2729        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2730        // Get the views  
2731        DataArrayView view_1 = tmp_1->getDefaultValue();        //Get the pointers to the actual data
2732        DataArrayView view_2 = tmp_2->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2733        // Get the pointers to the actual data        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2734        double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2735        // Compute a result for the default        // Compute a result for the default
2736        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2737        // Compute a result for each tag        // Compute a result for each tag
2738        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2739        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2740        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2741          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2742          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2743          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);  
2744          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2745        }        }
   
2746      }      }
2747      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2748    
# Line 2309  C_TensorBinaryOperation(Data const &arg_ Line 2755  C_TensorBinaryOperation(Data const &arg_
2755        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2756        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2757        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2758          res.requireWrite();
2759        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2760        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2761          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2762            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2763            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2764            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2765            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2766            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2767            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2768          }          }
2769        }        }
# Line 2337  C_TensorBinaryOperation(Data const &arg_ Line 2784  C_TensorBinaryOperation(Data const &arg_
2784    
2785        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2786        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2787        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();  
2788        // Get the pointers to the actual data        // Get the pointers to the actual data
2789        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2790        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2791        // Compute a result for the default        // Compute a result for the default
2792        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2793        // Compute a result for each tag        // Compute a result for each tag
2794        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2795        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2796        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2797          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2798          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2799          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);  
2800          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2801        }        }
2802    
# Line 2372  C_TensorBinaryOperation(Data const &arg_ Line 2814  C_TensorBinaryOperation(Data const &arg_
2814        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2815        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2816    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2817        // Get the pointers to the actual data        // Get the pointers to the actual data
2818        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2819        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2820        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2821    
2822        // Compute a result for the default        // Compute a result for the default
2823        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2824        // Merge the tags        // Merge the tags
# Line 2387  C_TensorBinaryOperation(Data const &arg_ Line 2826  C_TensorBinaryOperation(Data const &arg_
2826        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2827        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2828        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2829          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
2830        }        }
2831        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2832          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2833        }        }
2834        // Compute a result for each tag        // Compute a result for each tag
2835        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2836        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2837          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2838          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2839          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);  
2840          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2841        }        }
2842    
# Line 2416  C_TensorBinaryOperation(Data const &arg_ Line 2852  C_TensorBinaryOperation(Data const &arg_
2852        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2853        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2854        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2855          res.requireWrite();
2856        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2857        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2858          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2859          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2860          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2861            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2862            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2863            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2864            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2865            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2866          }          }
2867        }        }
2868    
2869      }      }
2870      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2871        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2872        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2873        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 2877  C_TensorBinaryOperation(Data const &arg_
2877        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2878        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2879        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2880          res.requireWrite();
2881        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2882        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2883          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2884            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2885            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2886            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2887            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2888            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2889            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2890          }          }
2891        }        }
# Line 2466  C_TensorBinaryOperation(Data const &arg_ Line 2903  C_TensorBinaryOperation(Data const &arg_
2903        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2904        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2905        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2906          res.requireWrite();
2907        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2908        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2909          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2910          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2911          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2912            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2913            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2914            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2915            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2916            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2917          }          }
2918        }        }
# Line 2491  C_TensorBinaryOperation(Data const &arg_ Line 2929  C_TensorBinaryOperation(Data const &arg_
2929        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2930        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2931        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2932          res.requireWrite();
2933        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2934        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2935          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2936            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2937            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2938            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2939            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2940            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2941            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2942            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2943          }          }
2944        }        }
# Line 2521  Data Line 2960  Data
2960  C_TensorUnaryOperation(Data const &arg_0,  C_TensorUnaryOperation(Data const &arg_0,
2961                         UnaryFunction operation)                         UnaryFunction operation)
2962  {  {
2963      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2964      {
2965         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2966      }
2967      if (arg_0.isLazy())
2968      {
2969         throw DataException("Error - Operations not permitted on lazy data.");
2970      }
2971    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2972    Data arg_0_Z = Data(arg_0);    Data arg_0_Z = Data(arg_0);
2973    
2974    // Get rank and shape of inputs    // Get rank and shape of inputs
2975    int rank0 = arg_0_Z.getDataPointRank();    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
   DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();  
2976    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2977    
2978    // Declare output Data object    // Declare output Data object
# Line 2534  C_TensorUnaryOperation(Data const &arg_0 Line 2980  C_TensorUnaryOperation(Data const &arg_0
2980    
2981    if (arg_0_Z.isConstant()) {    if (arg_0_Z.isConstant()) {
2982      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2983      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2984      double *ptr_2 = &((res.getPointDataView().getData())[0]);      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2985      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2986    }    }
2987    else if (arg_0_Z.isTagged()) {    else if (arg_0_Z.isTagged()) {
# Line 2548  C_TensorUnaryOperation(Data const &arg_0 Line 2994  C_TensorUnaryOperation(Data const &arg_0
2994      res.tag();      res.tag();
2995      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2996    
     // Get the views  
     DataArrayView view_0 = tmp_0->getDefaultValue();  
     DataArrayView view_2 = tmp_2->getDefaultValue();  
2997      // Get the pointers to the actual data      // Get the pointers to the actual data
2998      double *ptr_0 = &((view_0.getData())[0]);      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2999      double *ptr_2 = &((view_2.getData())[0]);      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
3000      // Compute a result for the default      // Compute a result for the default
3001      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3002      // Compute a result for each tag      // Compute a result for each tag
3003      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3004      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3005      for (i=lookup_0.begin();i!=lookup_0.end();i++) {      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3006        tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());        tmp_2->addTag(i->first);
3007        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
3008        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);  
3009        tensor_unary_operation(size0, ptr_0, ptr_2, operation);        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3010      }      }
3011    
# Line 2580  C_TensorUnaryOperation(Data const &arg_0 Line 3021  C_TensorUnaryOperation(Data const &arg_0
3021      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3022      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3023      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3024        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {      dataPointNo_0=0;
3025    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3026          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3027          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3028          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3029          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3030          tensor_unary_operation(size0, ptr_0, ptr_2, operation);          tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3031        }  //      }
3032      }      }
   
3033    }    }
3034    else {    else {
3035      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.2742

  ViewVC Help
Powered by ViewVC 1.1.26