/[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 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2010 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    
31  extern "C" {  extern "C" {
32  #include "DataC.h"  #include "DataC.h"
33  /* #include "paso/Paso.h" doesn't belong in this file...causes trouble for BruceFactory.cpp */  //#include <omp.h>
34  }  }
35    
36    #ifdef _OPENMP
37    #include <omp.h>
38    #endif
39    
40  #include "esysmpi.h"  #include "esysmpi.h"
41  #include <string>  #include <string>
42  #include <algorithm>  #include <algorithm>
43    #include <sstream>
44    
45  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
46  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
47  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
 #include <boost/python/numeric.hpp>  
48    
49  namespace escript {  namespace escript {
50    
# Line 47  namespace escript { Line 53  namespace escript {
53  class DataConstant;  class DataConstant;
54  class DataTagged;  class DataTagged;
55  class DataExpanded;  class DataExpanded;
56    class DataLazy;
57    
58  /**  /**
59     \brief     \brief
60     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
61    
62     Description:     Description:
63     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
64     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
65     the object created for the lifetime of the object.     of the Data object.
66     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
67     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
68       Doing so will lead to invalid memory access.
69       This should not affect any methods exposed via boost::python.
70  */  */
71  class Data {  class Data {
72    
# Line 101  class Data { Line 109  class Data {
109         const FunctionSpace& what);         const FunctionSpace& what);
110    
111    /**    /**
112       \brief      \brief Copy Data from an existing vector
113       Constructor which copies data from a DataArrayView.    */
114    
      \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.  
   */  
115    ESCRIPT_DLL_API    ESCRIPT_DLL_API
116    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
117         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
118         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
119                     bool expanded=false);
120    
121    /**    /**
122       \brief       \brief
123       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
124    
125       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
126       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 128  class Data { Line 131  class Data {
131    */    */
132    ESCRIPT_DLL_API    ESCRIPT_DLL_API
133    Data(double value,    Data(double value,
134         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
135         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
136         bool expanded=false);         bool expanded=false);
137    
# Line 141  class Data { Line 144  class Data {
144    */    */
145    ESCRIPT_DLL_API    ESCRIPT_DLL_API
146    Data(const Data& inData,    Data(const Data& inData,
147         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);  
148    
149    /**    /**
150       \brief       \brief
151       Constructor which copies data from a python numarray.       Constructor which copies data from any object that can be treated like a python array/sequence.
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   ESCRIPT_DLL_API  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from any object that can be converted into  
      a python numarray.  
152    
153       \param value - Input - Input data.       \param value - Input - Input data.
154       \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 164  class Data {
164    /**    /**
165       \brief       \brief
166       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
167       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
168       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
169    
170       \param value - Input - Input data.       \param value - Input - Input data.
171       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
# Line 217  class Data { Line 183  class Data {
183         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
184         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
185         bool expanded=false);         bool expanded=false);
186    
187    
188    
189      /**
190        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
191        Once you have passed the pointer, do not delete it.
192      */
193      ESCRIPT_DLL_API
194      explicit Data(DataAbstract* underlyingdata);
195    
196      /**
197        \brief Create a Data based on the supplied DataAbstract
198      */
199      ESCRIPT_DLL_API
200      explicit Data(DataAbstract_ptr underlyingdata);
201    
202    /**    /**
203       \brief       \brief
204       Destructor       Destructor
# Line 225  class Data { Line 207  class Data {
207    ~Data();    ~Data();
208    
209    /**    /**
210       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
211    */    */
212    ESCRIPT_DLL_API    ESCRIPT_DLL_API
213    void    void
214    copy(const Data& other);    copy(const Data& other);
215    
216    /**    /**
217         \brief Return a pointer to a deep copy of this object.
218      */
219      ESCRIPT_DLL_API
220      Data
221      copySelf();
222    
223    
224      /**
225         \brief produce a delayed evaluation version of this Data.
226      */
227      ESCRIPT_DLL_API
228      Data
229      delay();
230    
231      /**
232         \brief convert the current data into lazy data.
233      */
234      ESCRIPT_DLL_API
235      void
236      delaySelf();
237    
238    
239      /**
240       Member access methods.       Member access methods.
241    */    */
242    
# Line 247  class Data { Line 251  class Data {
251    
252    /**    /**
253       \brief       \brief
254       Returns trueif the data object is protected against update       Returns true, if the data object is protected against update
255    
256    */    */
257    ESCRIPT_DLL_API    ESCRIPT_DLL_API
258    bool    bool
259    isProtected() const;    isProtected() const;
260    
261    /**  
262       \brief  /**
263       Return the values of a data point on this process     \brief
264    */     Return the value of a data point as a python tuple.
265    */
266    ESCRIPT_DLL_API    ESCRIPT_DLL_API
267    const boost::python::numeric::array    const boost::python::object
268    getValueOfDataPoint(int dataPointNo);    getValueOfDataPointAsTuple(int dataPointNo);
269    
270    /**    /**
271       \brief       \brief
# Line 272  class Data { Line 277  class Data {
277    
278    /**    /**
279       \brief       \brief
280       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
281    */    */
282    ESCRIPT_DLL_API    ESCRIPT_DLL_API
283    void    void
284    setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array&);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
285    
286    /**    /**
287       \brief       \brief
# Line 287  class Data { Line 292  class Data {
292    setValueOfDataPoint(int dataPointNo, const double);    setValueOfDataPoint(int dataPointNo, const double);
293    
294    /**    /**
295       \brief       \brief Return a data point across all processors as a python tuple.
      Return the value of the specified data-point across all processors  
296    */    */
297    ESCRIPT_DLL_API    ESCRIPT_DLL_API
298    const boost::python::numeric::array    const boost::python::object
299    getValueOfGlobalDataPoint(int procNo, int dataPointNo);    getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
300    
301    /**    /**
302       \brief       \brief
303       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
304    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
305    */    */
306    ESCRIPT_DLL_API    ESCRIPT_DLL_API
307    int    int
# Line 313  class Data { Line 315  class Data {
315    escriptDataC    escriptDataC
316    getDataC();    getDataC();
317    
318    
319    
320    /**    /**
321       \brief       \brief
322       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 321  class Data { Line 325  class Data {
325    escriptDataC    escriptDataC
326    getDataC() const;    getDataC() const;
327    
   /**  
      \brief  
      Write the data as a string.  
   */  
   ESCRIPT_DLL_API  
   inline  
   std::string  
   toString() const  
   {  
     return m_data->toString();  
   }  
328    
329    /**    /**
330       \brief       \brief
331       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.  
332    */    */
333    ESCRIPT_DLL_API    ESCRIPT_DLL_API
334    inline    std::string
335    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
336    
337    /**    /**
338       \brief       \brief
# Line 360  class Data { Line 347  class Data {
347       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
348       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
349       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
350    */    */
351    ESCRIPT_DLL_API    ESCRIPT_DLL_API
352    void    void
353    tag();    tag();
354    
355    /**    /**
356        \brief If this data is lazy, then convert it to ready data.
357        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
358      */
359      ESCRIPT_DLL_API
360      void
361      resolve();
362    
363    
364      /**
365       \brief Ensures data is ready for write access.
366      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
367      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
368      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
369      Doing so might introduce additional sharing.
370      */
371      ESCRIPT_DLL_API
372      void
373      requireWrite();
374    
375      /**
376       \brief       \brief
377       Return true if this Data is expanded.       Return true if this Data is expanded.
378         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
379    */    */
380    ESCRIPT_DLL_API    ESCRIPT_DLL_API
381    bool    bool
# Line 376  class Data { Line 383  class Data {
383    
384    /**    /**
385       \brief       \brief
386         Return true if this Data is expanded or resolves to expanded.
387         That is, if it has a separate value for each datapoint in the sample.
388      */
389      ESCRIPT_DLL_API
390      bool
391      actsExpanded() const;
392      
393    
394      /**
395         \brief
396       Return true if this Data is tagged.       Return true if this Data is tagged.
397    */    */
398    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 391  class Data { Line 408  class Data {
408    isConstant() const;    isConstant() const;
409    
410    /**    /**
411         \brief Return true if this Data is lazy.
412      */
413      ESCRIPT_DLL_API
414      bool
415      isLazy() const;
416    
417      /**
418         \brief Return true if this data is ready.
419      */
420      ESCRIPT_DLL_API
421      bool
422      isReady() const;
423    
424      /**
425       \brief       \brief
426       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
427    contains datapoints.
428    */    */
429    ESCRIPT_DLL_API    ESCRIPT_DLL_API
430    bool    bool
# Line 424  class Data { Line 456  class Data {
456    */    */
457    ESCRIPT_DLL_API    ESCRIPT_DLL_API
458    inline    inline
459    const AbstractDomain&  //   const AbstractDomain&
460      const_Domain_ptr
461    getDomain() const    getDomain() const
462    {    {
463       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
464    }    }
465    
466    
467      /**
468         \brief
469         Return the domain.
470         TODO: For internal use only.   This should be removed.
471      */
472      ESCRIPT_DLL_API
473      inline
474    //   const AbstractDomain&
475      Domain_ptr
476      getDomainPython() const
477      {
478         return getFunctionSpace().getDomainPython();
479      }
480    
481    /**    /**
482       \brief       \brief
483       Return a copy of the domain.       Return a copy of the domain.
# Line 444  class Data { Line 492  class Data {
492    */    */
493    ESCRIPT_DLL_API    ESCRIPT_DLL_API
494    inline    inline
495    int    unsigned int
496    getDataPointRank() const    getDataPointRank() const
497    {    {
498      return m_data->getPointDataView().getRank();      return m_data->getRank();
499    }    }
500    
501    /**    /**
# Line 484  class Data { Line 532  class Data {
532    {    {
533      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
534    }    }
535    
536    
537      /**
538        \brief
539        Return the number of values in the shape for this object.
540      */
541      ESCRIPT_DLL_API
542      int
543      getNoValues() const
544      {
545        return m_data->getNoValues();
546      }
547    
548    
549    /**    /**
550       \brief       \brief
551       dumps the object into a netCDF file       dumps the object into a netCDF file
# Line 491  class Data { Line 553  class Data {
553    ESCRIPT_DLL_API    ESCRIPT_DLL_API
554    void    void
555    dump(const std::string fileName) const;    dump(const std::string fileName) const;
556    
557     /**
558      \brief returns the values of the object as a list of tuples (one for each datapoint).
559    
560      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
561    If false, the result is a list of scalars [1, 2, ...]
562     */
563      ESCRIPT_DLL_API
564      const boost::python::object
565      toListOfTuples(bool scalarastuple=true);
566    
567    
568     /**
569        \brief
570        Return the sample data for the given sample no. This is not the
571        preferred interface but is provided for use by C code.
572        \param sampleNo - Input - the given sample no.
573        \return pointer to the sample data.
574    */
575      ESCRIPT_DLL_API
576      inline
577      const DataAbstract::ValueType::value_type*
578      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo);
579    
580    
581    /**    /**
582       \brief       \brief
583       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
584       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
585       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
586         \return pointer to the sample data.
587    */    */
588    ESCRIPT_DLL_API    ESCRIPT_DLL_API
589    inline    inline
590    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
591    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
592    {  
     return m_data->getSampleData(sampleNo);  
   }  
593    
594    /**    /**
595       \brief       \brief
# Line 521  class Data { Line 607  class Data {
607    
608    /**    /**
609       \brief       \brief
610       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.  
611       \param sampleNo - Input -       \param sampleNo - Input -
612       \param dataPointNo - Input -       \param dataPointNo - Input -
613    */    */
614    ESCRIPT_DLL_API    ESCRIPT_DLL_API
615      DataTypes::ValueType::const_reference
616      getDataPointRO(int sampleNo, int dataPointNo);
617    
618      /**
619         \brief
620         Return a reference into the DataVector which points to the specified data point.
621         \param sampleNo - Input -
622         \param dataPointNo - Input -
623      */
624      ESCRIPT_DLL_API
625      DataTypes::ValueType::reference
626      getDataPointRW(int sampleNo, int dataPointNo);
627    
628    
629    
630      /**
631         \brief
632         Return the offset for the given sample and point within the sample
633      */
634      ESCRIPT_DLL_API
635    inline    inline
636    DataArrayView    DataTypes::ValueType::size_type
637    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
638                 int dataPointNo)                 int dataPointNo)
639    {    {
640                  return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
641    }    }
642    
643    /**    /**
# Line 541  class Data { Line 645  class Data {
645       Return a reference to the data point shape.       Return a reference to the data point shape.
646    */    */
647    ESCRIPT_DLL_API    ESCRIPT_DLL_API
648    const DataArrayView::ShapeType&    inline
649    getDataPointShape() const;    const DataTypes::ShapeType&
650      getDataPointShape() const
651      {
652        return m_data->getShape();
653      }
654    
655    /**    /**
656       \brief       \brief
# Line 566  class Data { Line 674  class Data {
674       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
675    */    */
676    ESCRIPT_DLL_API    ESCRIPT_DLL_API
677    DataArrayView::ValueType::size_type    DataTypes::ValueType::size_type
678    getLength() const;    getLength() const;
679    
680      /**
681      \brief Return true if this object contains no samples.
682      This is not the same as isEmpty()
683      */
684      ESCRIPT_DLL_API
685      bool
686      hasNoSamples() const
687      {
688        return getLength()==0;
689      }
690    
691    /**    /**
692       \brief       \brief
693       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
694       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
695       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.
696       \param tagKey - Input - Integer key.       \param name - Input - name of tag.
697       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
     ==>*  
698    */    */
699    ESCRIPT_DLL_API    ESCRIPT_DLL_API
700    void    void
# Line 605  class Data { Line 721  class Data {
721       object to type DataTagged if it is constant.       object to type DataTagged if it is constant.
722    
723       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
724         \param pointshape - Input - The shape of the value parameter
725       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
726      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
727    */    */
728    ESCRIPT_DLL_API    ESCRIPT_DLL_API
729    void    void
730    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
731                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
732                            const DataTypes::ValueType& value,
733                int dataOffset=0);
734    
735    
736    
737    /**    /**
738      \brief      \brief
# Line 644  class Data { Line 765  class Data {
765    ESCRIPT_DLL_API    ESCRIPT_DLL_API
766    Data    Data
767    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
768    
769    
770      ESCRIPT_DLL_API
771      Data
772      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
773                           double undef, Data& B, double Bmin, double Bstep,bool check_boundaries);
774    
775      ESCRIPT_DLL_API
776      Data
777      interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
778                           double undef,bool check_boundaries);
779    
780    
781    
782    
783      ESCRIPT_DLL_API
784      Data
785      interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
786                            Data& B, double Bmin, double Bstep, double undef,bool check_boundaries);
787    
788      ESCRIPT_DLL_API
789      Data
790      interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
791                            double undef,bool check_boundaries);
792    
793    /**    /**
794       \brief       \brief
795       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 805  class Data {
805    grad() const;    grad() const;
806    
807    /**    /**
808       \brief      \brief
809       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
      *  
810    */    */
811    ESCRIPT_DLL_API    ESCRIPT_DLL_API
812    boost::python::numeric::array    boost::python::object
813    integrate() const;    integrateToTuple_const() const;
814    
815    
816      /**
817        \brief
818         Calculate the integral over the function space domain as a python tuple.
819      */
820      ESCRIPT_DLL_API
821      boost::python::object
822      integrateToTuple();
823    
824    
825    
826    /**    /**
827       \brief       \brief
# Line 732  class Data { Line 888  class Data {
888    /**    /**
889       \brief       \brief
890       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
891       *  
892         The method is not const because lazy data needs to be expanded before Lsup can be computed.
893         The _const form can be used when the Data object is const, however this will only work for
894         Data which is not Lazy.
895    
896         For Data which contain no samples (or tagged Data for which no tags in use have a value)
897         zero is returned.
898    */    */
899    ESCRIPT_DLL_API    ESCRIPT_DLL_API
900    double    double
901    Lsup() const;    Lsup();
902    
903      ESCRIPT_DLL_API
904      double
905      Lsup_const() const;
906    
907    
908    /**    /**
909       \brief       \brief
910       Return the maximum value of this Data object.       Return the maximum value of this Data object.
911       *  
912         The method is not const because lazy data needs to be expanded before sup can be computed.
913         The _const form can be used when the Data object is const, however this will only work for
914         Data which is not Lazy.
915    
916         For Data which contain no samples (or tagged Data for which no tags in use have a value)
917         a large negative value is returned.
918    */    */
919    ESCRIPT_DLL_API    ESCRIPT_DLL_API
920    double    double
921    sup() const;    sup();
922    
923      ESCRIPT_DLL_API
924      double
925      sup_const() const;
926    
927    
928    /**    /**
929       \brief       \brief
930       Return the minimum value of this Data object.       Return the minimum value of this Data object.
931       *  
932         The method is not const because lazy data needs to be expanded before inf can be computed.
933         The _const form can be used when the Data object is const, however this will only work for
934         Data which is not Lazy.
935    
936         For Data which contain no samples (or tagged Data for which no tags in use have a value)
937         a large positive value is returned.
938    */    */
939    ESCRIPT_DLL_API    ESCRIPT_DLL_API
940    double    double
941    inf() const;    inf();
942    
943      ESCRIPT_DLL_API
944      double
945      inf_const() const;
946    
947    
948    
949    /**    /**
950       \brief       \brief
# Line 786  class Data { Line 976  class Data {
976    /**    /**
977       \brief       \brief
978       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
979       the minimum value in this Data object.       the minimum component value in this Data object.
980         \note If you are working in python, please consider using Locator
981    instead of manually manipulating process and point IDs.
982    */    */
983    ESCRIPT_DLL_API    ESCRIPT_DLL_API
984    const boost::python::tuple    const boost::python::tuple
985    minGlobalDataPoint() const;    minGlobalDataPoint() const;
986    
987      /**
988         \brief
989         Return the (sample number, data-point number) of the data point with
990         the minimum component value in this Data object.
991         \note If you are working in python, please consider using Locator
992    instead of manually manipulating process and point IDs.
993      */
994    ESCRIPT_DLL_API    ESCRIPT_DLL_API
995    void    const boost::python::tuple
996    calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;    maxGlobalDataPoint() const;
997    
998    
999    
1000    /**    /**
1001       \brief       \brief
1002       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 1297  class Data {
1297    void    void
1298    saveVTK(std::string fileName) const;    saveVTK(std::string fileName) const;
1299    
1300    
1301    
1302    /**    /**
1303       \brief       \brief
1304       Overloaded operator +=       Overloaded operator +=
# Line 1143  class Data { Line 1347  class Data {
1347    Data& operator/=(const boost::python::object& right);    Data& operator/=(const boost::python::object& right);
1348    
1349    /**    /**
1350        \brief return inverse of matricies.
1351      */
1352      ESCRIPT_DLL_API
1353      Data
1354      matrixInverse() const;
1355    
1356      /**
1357       \brief       \brief
1358       Returns true if this can be interpolated to functionspace.       Returns true if this can be interpolated to functionspace.
1359    */    */
# Line 1212  class Data { Line 1423  class Data {
1423    */    */
1424    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1425    Data    Data
1426    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1427    
1428    /**    /**
1429       \brief       \brief
# Line 1225  class Data { Line 1436  class Data {
1436    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1437    void    void
1438    setSlice(const Data& value,    setSlice(const Data& value,
1439             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Archive the current Data object to the given file.  
      \param fileName - Input - file to archive to.  
   */  
   ESCRIPT_DLL_API  
   void  
   archiveData(const std::string fileName);  
   
   /**  
      \brief  
      Extract the Data object archived in the given file, overwriting  
      the current Data object.  
      Note - the current object must be of type DataEmpty.  
      \param fileName - Input - file to extract from.  
      \param fspace - Input - a suitable FunctionSpace descibing the data.  
   */  
   ESCRIPT_DLL_API  
   void  
   extractData(const std::string fileName,  
               const FunctionSpace& fspace);  
   
1440    
1441    /**    /**
1442       \brief       \brief
# Line 1290  class Data { Line 1478  class Data {
1478    /**    /**
1479       \brief       \brief
1480       return the object produced by the factory, which is a DataConstant or DataExpanded       return the object produced by the factory, which is a DataConstant or DataExpanded
1481        TODO Ownership of this object should be explained in doco.
1482    */    */
1483    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1484          DataAbstract*          DataAbstract*
1485          borrowData(void) const;          borrowData(void) const;
1486    
1487      ESCRIPT_DLL_API
1488            DataAbstract_ptr
1489            borrowDataPtr(void) const;
1490    
1491      ESCRIPT_DLL_API
1492            DataReady_ptr
1493            borrowReadyPtr(void) const;
1494    
1495    
1496    
1497      /**
1498         \brief
1499         Return a pointer to the beginning of the datapoint at the specified offset.
1500         TODO Eventually these should be inlined.
1501         \param i - position(offset) in the underlying datastructure
1502      */
1503    
1504      ESCRIPT_DLL_API
1505            DataTypes::ValueType::const_reference
1506            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1507    
1508    
1509      ESCRIPT_DLL_API
1510            DataTypes::ValueType::reference
1511            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1512    
1513    
1514    
1515   protected:   protected:
1516    
1517   private:   private:
1518    
1519    template <class BinaryOp>
1520      double
1521    #ifdef PASO_MPI
1522      lazyAlgWorker(double init, MPI_Op mpiop_type);
1523    #else
1524      lazyAlgWorker(double init);
1525    #endif
1526    
1527      double
1528      LsupWorker() const;
1529    
1530      double
1531      supWorker() const;
1532    
1533      double
1534      infWorker() const;
1535    
1536      boost::python::object
1537      integrateWorker() const;
1538    
1539      void
1540      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1541    
1542      void
1543      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1544    
1545      // For internal use in Data.cpp only!
1546      // other uses should call the main entry points and allow laziness
1547      Data
1548      minval_nonlazy() const;
1549    
1550      // For internal use in Data.cpp only!
1551      Data
1552      maxval_nonlazy() const;
1553    
1554    
1555    /**    /**
1556       \brief       \brief
1557       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1370  class Data { Line 1623  class Data {
1623       \brief       \brief
1624       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1625    */    */
1626    template <class IValueType>  
1627    void    void
1628    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1629             const DataTypes::ShapeType& shape,
1630               const FunctionSpace& what,               const FunctionSpace& what,
1631               bool expanded);               bool expanded);
1632    
1633      void
1634      initialise(const WrappedArray& value,
1635                     const FunctionSpace& what,
1636                     bool expanded);
1637    
1638    //    //
1639    // flag to protect the data object against any update    // flag to protect the data object against any update
1640    bool m_protected;    bool m_protected;
1641      mutable bool m_shared;
1642      bool m_lazy;
1643    
1644    //    //
1645    // pointer to the actual data object    // pointer to the actual data object
1646    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1647      DataAbstract_ptr m_data;
1648    
1649  };  // If possible please use getReadyPtr instead.
1650    // But see warning below.
1651      const DataReady*
1652      getReady() const;
1653    
1654  template <class IValueType>    DataReady*
1655  void    getReady();
1656  Data::initialise(const IValueType& value,  
1657                   const FunctionSpace& what,  
1658                   bool expanded)  // Be wary of using this for local operations since it (temporarily) increases reference count.
1659  {  // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1660    // getReady() instead
1661      DataReady_ptr
1662      getReadyPtr();
1663    
1664      const_DataReady_ptr
1665      getReadyPtr() const;
1666    
1667    
1668      /**
1669       \brief Update the Data's shared flag
1670       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1671       For internal use only.
1672      */
1673      void updateShareStatus(bool nowshared) const
1674      {
1675        m_shared=nowshared;     // m_shared is mutable
1676      }
1677    
1678      // In the isShared() method below:
1679      // A problem would occur if m_data (the address pointed to) were being modified
1680      // while the call m_data->is_shared is being executed.
1681      //
1682      // Q: So why do I think this code can be thread safe/correct?
1683      // A: We need to make some assumptions.
1684      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1685      // 2. We assume that no constructions or assignments which will share previously unshared
1686      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1687    //    //
1688    // 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
1689    // 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.
1690    // 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.
1691    // within the shared_ptr constructor.    // For any threads executing before the flag switches they will assume the object is still shared.
1692    if (expanded) {    bool isShared() const
1693      DataAbstract* temp=new DataExpanded(value,what);    {
1694      boost::shared_ptr<DataAbstract> temp_data(temp);      return m_shared;
1695      m_data=temp_data;  /*  if (m_shared) return true;
1696    } else {      if (m_data->isShared())        
1697      DataAbstract* temp=new DataConstant(value,what);      {                  
1698      boost::shared_ptr<DataAbstract> temp_data(temp);          updateShareStatus(true);
1699      m_data=temp_data;          return true;
1700        }
1701        return false;*/
1702      }
1703    
1704      void forceResolve()
1705      {
1706        if (isLazy())
1707        {
1708            #ifdef _OPENMP
1709            if (omp_in_parallel())
1710            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1711            throw DataException("Please do not call forceResolve() in a parallel region.");
1712            }
1713            #endif
1714            resolve();
1715        }
1716      }
1717    
1718      /**
1719      \brief if another object is sharing out member data make a copy to work with instead.
1720      This code should only be called from single threaded sections of code.
1721      */
1722      void exclusiveWrite()
1723      {
1724    #ifdef _OPENMP
1725      if (omp_in_parallel())
1726      {
1727    // *((int*)0)=17;
1728        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1729      }
1730    #endif
1731        forceResolve();
1732        if (isShared())
1733        {
1734            DataAbstract* t=m_data->deepCopy();
1735            set_m_data(DataAbstract_ptr(t));
1736        }
1737    }    }
1738    
1739      /**
1740      \brief checks if caller can have exclusive write to the object
1741      */
1742      void checkExclusiveWrite()
1743      {
1744        if  (isLazy() || isShared())
1745        {
1746            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1747        }
1748      }
1749    
1750      /**
1751      \brief Modify the data abstract hosted by this Data object
1752      For internal use only.
1753      Passing a pointer to null is permitted (do this in the destructor)
1754      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1755      */
1756      void set_m_data(DataAbstract_ptr p);
1757    
1758      friend class DataAbstract;        // To allow calls to updateShareStatus
1759    #ifdef IKNOWWHATIMDOING
1760      friend Data applyBinaryCFunction(boost::python::object cfunc, boost::python::tuple shape, escript::Data& d, escript::Data& e);
1761    #endif
1762    };
1763    
1764    
1765    #ifdef IKNOWWHATIMDOING
1766    ESCRIPT_DLL_API
1767    Data
1768    applyBinaryCFunction(boost::python::object func, boost::python::tuple shape, escript::Data& d, escript::Data& e);
1769    #endif
1770    
1771    
1772    }   // end namespace escript
1773    
1774    
1775    // No, this is not supposed to be at the top of the file
1776    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1777    // so that I can dynamic cast between them below.
1778    #include "DataReady.h"
1779    #include "DataLazy.h"
1780    
1781    namespace escript
1782    {
1783    
1784    inline
1785    const DataReady*
1786    Data::getReady() const
1787    {
1788       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1789       EsysAssert((dr!=0), "Error - casting to DataReady.");
1790       return dr;
1791    }
1792    
1793    inline
1794    DataReady*
1795    Data::getReady()
1796    {
1797       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1798       EsysAssert((dr!=0), "Error - casting to DataReady.");
1799       return dr;
1800    }
1801    
1802    // Be wary of using this for local operations since it (temporarily) increases reference count.
1803    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1804    // getReady() instead
1805    inline
1806    DataReady_ptr
1807    Data::getReadyPtr()
1808    {
1809       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1810       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1811       return dr;
1812  }  }
1813    
1814    
1815    inline
1816    const_DataReady_ptr
1817    Data::getReadyPtr() const
1818    {
1819       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1820       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1821       return dr;
1822    }
1823    
1824    inline
1825    DataAbstract::ValueType::value_type*
1826    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1827    {
1828       if (isLazy())
1829       {
1830        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1831       }
1832       return getReady()->getSampleDataRW(sampleNo);
1833    }
1834    
1835    inline
1836    const DataAbstract::ValueType::value_type*
1837    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo)
1838    {
1839       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1840       if (l!=0)
1841       {
1842        size_t offset=0;
1843        const DataTypes::ValueType* res=l->resolveSample(sampleNo,offset);
1844        return &((*res)[offset]);
1845       }
1846       return getReady()->getSampleDataRO(sampleNo);
1847    }
1848    
1849    
1850    
1851    /**
1852       Modify a filename for MPI parallel output to multiple files
1853    */
1854    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1855    
1856  /**  /**
1857     Binary Data object operators.     Binary Data object operators.
1858  */  */
# Line 1519  ESCRIPT_DLL_API std::ostream& operator<< Line 1964  ESCRIPT_DLL_API std::ostream& operator<<
1964  /**  /**
1965    \brief    \brief
1966    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1967    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1968    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1969    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1970    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1971  */  */
1972  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1973  Data  Data
1974  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1975                       Data& arg1,                       Data& arg_1,
1976                       int axis_offset=0,                       int axis_offset=0,
1977                       int transpose=0);                       int transpose=0);
1978    
   
 /**  
   \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);  
   
1979  /**  /**
1980    \brief    \brief
1981    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 1989  Data::binaryOp(const Data& right,
1989  {  {
1990     //     //
1991     // 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
1992     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1993       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.");
1994     }     }
1995    
1996       if (isLazy() || right.isLazy())
1997       {
1998         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1999       }
2000     //     //
2001     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
2002     Data tempRight(right);     Data tempRight(right);
2003    
2004     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
2005       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
2006         //         //
# Line 1582  Data::binaryOp(const Data& right, Line 2010  Data::binaryOp(const Data& right,
2010         //         //
2011         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2012         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2013         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2014           set_m_data(tempLeft.m_data);
2015       }       }
2016     }     }
2017     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1598  Data::binaryOp(const Data& right, Line 2027  Data::binaryOp(const Data& right,
2027       // of any data type       // of any data type
2028       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2029       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2030       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2031     } else if (isTagged()) {     } else if (isTagged()) {
2032       //       //
2033       // 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 2075  Data::algorithm(BinaryFunction operation
2075      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2076      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2077      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2078      } else if (isEmpty()) {
2079        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2080      } else if (isLazy()) {
2081        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2082      } else {
2083        throw DataException("Error - Data encapsulates an unknown type.");
2084    }    }
   return 0;  
2085  }  }
2086    
2087  /**  /**
# Line 1663  inline Line 2097  inline
2097  Data  Data
2098  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2099  {  {
2100    if (isExpanded()) {    if (isEmpty()) {
2101      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2102      }
2103      else if (isExpanded()) {
2104        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2105      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2106      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2107      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2108      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2109      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2110      return result;      return result;
2111    } else if (isTagged()) {    }
2112      else if (isTagged()) {
2113      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());  
2114      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2115      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2116        defval[0]=0;
2117        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2118      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2119      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2120    } else if (isConstant()) {    }
2121      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2122        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2123      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2124      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2125      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2126      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2127      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2128      return result;      return result;
2129      } else if (isLazy()) {
2130        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2131      } else {
2132        throw DataException("Error - Data encapsulates an unknown type.");
2133    }    }
   Data falseRetVal; // to keep compiler quiet  
   return falseRetVal;  
2134  }  }
2135    
2136    /**
2137      \brief
2138      Compute a tensor operation with two Data objects
2139      \param arg_0 - Input - Data object
2140      \param arg_1 - Input - Data object
2141      \param operation - Input - Binary op functor
2142    */
2143  template <typename BinaryFunction>  template <typename BinaryFunction>
2144    inline
2145  Data  Data
2146  C_TensorBinaryOperation(Data const &arg_0,  C_TensorBinaryOperation(Data const &arg_0,
2147                          Data const &arg_1,                          Data const &arg_1,
2148                          BinaryFunction operation)                          BinaryFunction operation)
2149  {  {
2150      if (arg_0.isEmpty() || arg_1.isEmpty())
2151      {
2152         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2153      }
2154      if (arg_0.isLazy() || arg_1.isLazy())
2155      {
2156         throw DataException("Error - Operations not permitted on lazy data.");
2157      }
2158    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2159    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
2160    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
# Line 1730  C_TensorBinaryOperation(Data const &arg_ Line 2176  C_TensorBinaryOperation(Data const &arg_
2176    // Get rank and shape of inputs    // Get rank and shape of inputs
2177    int rank0 = arg_0_Z.getDataPointRank();    int rank0 = arg_0_Z.getDataPointRank();
2178    int rank1 = arg_1_Z.getDataPointRank();    int rank1 = arg_1_Z.getDataPointRank();
2179    DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();    DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2180    DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();    DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2181    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2182    int size1 = arg_1_Z.getDataPointSize();    int size1 = arg_1_Z.getDataPointSize();
   
2183    // Declare output Data object    // Declare output Data object
2184    Data res;    Data res;
2185    
2186    if (shape0 == shape1) {    if (shape0 == shape1) {
   
2187      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2188        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2189        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2190        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2191        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2192    
2193        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2194      }      }
2195      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 2207  C_TensorBinaryOperation(Data const &arg_
2207    
2208        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2209        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2210        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2211        // Get the views  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2212        // Get the pointers to the actual data        // Get the pointers to the actual data
2213        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2214        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2215    
2216        // Compute a result for the default        // Compute a result for the default
2217        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2218        // Compute a result for each tag        // Compute a result for each tag
2219        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2220        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
2221        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2222          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2223          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2224          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2225          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2226          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2227        }        }
2228    
2229      }      }
2230      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
   
2231        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2232        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2233        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 2237  C_TensorBinaryOperation(Data const &arg_
2237        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2238        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2239        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2240          res.requireWrite();
2241        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2242        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2243          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2244            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2245            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2246            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2247            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2248            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2249            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2250          }          }
2251        }        }
2252    
2253      }      }
2254      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
   
2255        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2256        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2257    
# Line 1823  C_TensorBinaryOperation(Data const &arg_ Line 2265  C_TensorBinaryOperation(Data const &arg_
2265    
2266        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2267        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2268        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);  
2269        // 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();  
2270        // Get the pointers to the actual data        // Get the pointers to the actual data
2271        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2272        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2273        // Compute a result for the default        // Compute a result for the default
2274        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2275        // Compute a result for each tag        // Compute a result for each tag
2276        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2277        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
2278        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2279          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2280          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2281          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);  
2282          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2283        }        }
2284    
2285      }      }
2286      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
   
2287        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2288        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2289    
# Line 1858  C_TensorBinaryOperation(Data const &arg_ Line 2295  C_TensorBinaryOperation(Data const &arg_
2295        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2296        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2297    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2298        // Get the pointers to the actual data        // Get the pointers to the actual data
2299        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2300        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2301        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2302    
2303        // Compute a result for the default        // Compute a result for the default
2304        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2305        // Merge the tags        // Merge the tags
# Line 1873  C_TensorBinaryOperation(Data const &arg_ Line 2307  C_TensorBinaryOperation(Data const &arg_
2307        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2308        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2309        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2310          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
2311        }        }
2312        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2313          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2314        }        }
2315        // Compute a result for each tag        // Compute a result for each tag
2316        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2317        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2318          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);  
2319          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2320          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2321          double *ptr_0 = &view_0.getData(0);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2322          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2323          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2324        }        }
2325    
2326      }      }
2327      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
   
2328        // 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
2329        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2330        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 2334  C_TensorBinaryOperation(Data const &arg_
2334        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2335        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2336        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2337          res.requireWrite();
2338        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2339        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2340          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
2341          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2342          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2343            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2344            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2345            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2346            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2347            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2348          }          }
2349        }        }
2350    
2351      }      }
2352      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2353        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2354        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2355        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 2359  C_TensorBinaryOperation(Data const &arg_
2359        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2360        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2361        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2362          res.requireWrite();
2363        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2364        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2365          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2366            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2367            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2368            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);  
2369            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2370            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2371              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2372    
2373    
2374            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2375          }          }
2376        }        }
2377    
2378      }      }
2379      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
   
2380        // 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
2381        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2382        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 2386  C_TensorBinaryOperation(Data const &arg_
2386        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2387        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2388        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2389          res.requireWrite();
2390        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2391        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2392          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2393          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2394          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2395            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2396            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2397            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2398            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2399            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2400          }          }
2401        }        }
2402    
2403      }      }
2404      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
   
2405        // 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
2406        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2407        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 2411  C_TensorBinaryOperation(Data const &arg_
2411        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2412        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2413        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2414          res.requireWrite();
2415        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2416        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2417          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        dataPointNo_0=0;
2418    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2419            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2420            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2421            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2422            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2423            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2424            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2425            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2426          }  //       }
2427        }        }
2428    
2429      }      }
# Line 1995  C_TensorBinaryOperation(Data const &arg_ Line 2432  C_TensorBinaryOperation(Data const &arg_
2432      }      }
2433    
2434    } else if (0 == rank0) {    } else if (0 == rank0) {
   
2435      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2436        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2437        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2438        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2439        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2440        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2441      }      }
2442      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 2454  C_TensorBinaryOperation(Data const &arg_
2454    
2455        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2456        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2457        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2458        // Get the views  
2459        DataArrayView view_1 = tmp_1->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2460        DataArrayView view_2 = tmp_2->getDefaultValue();        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2461        // Get the pointers to the actual data  
       double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2462        // Compute a result for the default        // Compute a result for the default
2463        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2464        // Compute a result for each tag        // Compute a result for each tag
2465        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2466        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
2467        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2468          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2469          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2470          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);  
2471          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2472        }        }
2473    
# Line 2047  C_TensorBinaryOperation(Data const &arg_ Line 2479  C_TensorBinaryOperation(Data const &arg_
2479        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2480        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2481    
2482        int sampleNo_1,dataPointNo_1;        int sampleNo_1;
2483        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2484        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2485        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2486        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        const double *ptr_src = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2487          double ptr_0 = ptr_src[0];
2488          int size = size1*numDataPointsPerSample_1;
2489          res.requireWrite();
2490          #pragma omp parallel for private(sampleNo_1) schedule(static)
2491        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2492          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {  //        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2493            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,0);
2494            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,0);
2495            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);  //          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2496            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2497            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2498            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
2499    
2500          }  //        }
2501        }        }
2502    
2503      }      }
# Line 2080  C_TensorBinaryOperation(Data const &arg_ Line 2516  C_TensorBinaryOperation(Data const &arg_
2516    
2517        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2518        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2519        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2520        // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2521        // Get the pointers to the actual data        // Get the pointers to the actual data
2522        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2523        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2524    
2525    
2526        // Compute a result for the default        // Compute a result for the default
2527        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2528        // Compute a result for each tag        // Compute a result for each tag
2529        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2530        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
2531        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2532          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2533          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2534          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2535          double *ptr_0 = &view_0.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2536          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2537        }        }
2538    
# Line 2115  C_TensorBinaryOperation(Data const &arg_ Line 2550  C_TensorBinaryOperation(Data const &arg_
2550        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2551        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2552    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2553        // Get the pointers to the actual data        // Get the pointers to the actual data
2554        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2555        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2556        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2557    
2558        // Compute a result for the default        // Compute a result for the default
2559        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2560        // Merge the tags        // Merge the tags
# Line 2130  C_TensorBinaryOperation(Data const &arg_ Line 2562  C_TensorBinaryOperation(Data const &arg_
2562        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2563        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2564        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2565          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
2566        }        }
2567        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2568          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2569        }        }
2570        // Compute a result for each tag        // Compute a result for each tag
2571        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2572        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2573          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2574          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2575          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2576          double *ptr_0 = &view_0.getData(0);  
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
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        }        }
2579    
# Line 2159  C_TensorBinaryOperation(Data const &arg_ Line 2589  C_TensorBinaryOperation(Data const &arg_
2589        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2590        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2591        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2592          res.requireWrite();
2593        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2594        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2595          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
2596          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2597          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2598            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2599            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2600            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2601            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2602            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2603          }          }
2604        }        }
2605    
2606      }      }
2607      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2608        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2609        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2610        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 2614  C_TensorBinaryOperation(Data const &arg_
2614        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2615        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2616        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2617          res.requireWrite();
2618        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2619        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2620          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2621            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2622            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2623            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2624            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2625            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2626            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2627          }          }
2628        }        }
# Line 2209  C_TensorBinaryOperation(Data const &arg_ Line 2640  C_TensorBinaryOperation(Data const &arg_
2640        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2641        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2642        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2643          res.requireWrite();
2644        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2645        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2646          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2647          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2648          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2649            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2650            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2651            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2652            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2653            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2654          }          }
2655        }        }
# Line 2234  C_TensorBinaryOperation(Data const &arg_ Line 2666  C_TensorBinaryOperation(Data const &arg_
2666        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2667        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2668        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2669          res.requireWrite();
2670        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2671        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2672          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2673            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2674            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2675            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2676            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2677            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2678            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2679            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2680          }          }
2681        }        }
# Line 2253  C_TensorBinaryOperation(Data const &arg_ Line 2686  C_TensorBinaryOperation(Data const &arg_
2686      }      }
2687    
2688    } else if (0 == rank1) {    } else if (0 == rank1) {
   
2689      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2690        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2691        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2692        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2693        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2694        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2695      }      }
2696      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 2708  C_TensorBinaryOperation(Data const &arg_
2708    
2709        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2710        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2711        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2712        // Get the views  
2713        DataArrayView view_1 = tmp_1->getDefaultValue();        //Get the pointers to the actual data
2714        DataArrayView view_2 = tmp_2->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2715        // Get the pointers to the actual data        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2716        double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2717        // Compute a result for the default        // Compute a result for the default
2718        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2719        // Compute a result for each tag        // Compute a result for each tag
2720        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2721        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
2722        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2723          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2724          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2725          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);  
2726          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2727        }        }
   
2728      }      }
2729      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2730    
# Line 2309  C_TensorBinaryOperation(Data const &arg_ Line 2737  C_TensorBinaryOperation(Data const &arg_
2737        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2738        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2739        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2740          res.requireWrite();
2741        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2742        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2743          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2744            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2745            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2746            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2747            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2748            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2749            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2750          }          }
2751        }        }
# Line 2337  C_TensorBinaryOperation(Data const &arg_ Line 2766  C_TensorBinaryOperation(Data const &arg_
2766    
2767        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2768        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2769        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();  
2770        // Get the pointers to the actual data        // Get the pointers to the actual data
2771        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2772        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2773        // Compute a result for the default        // Compute a result for the default
2774        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2775        // Compute a result for each tag        // Compute a result for each tag
2776        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2777        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
2778        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2779          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2780          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2781          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);  
2782          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2783        }        }
2784    
# Line 2372  C_TensorBinaryOperation(Data const &arg_ Line 2796  C_TensorBinaryOperation(Data const &arg_
2796        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2797        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2798    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2799        // Get the pointers to the actual data        // Get the pointers to the actual data
2800        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2801        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2802        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2803    
2804        // Compute a result for the default        // Compute a result for the default
2805        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2806        // Merge the tags        // Merge the tags
# Line 2387  C_TensorBinaryOperation(Data const &arg_ Line 2808  C_TensorBinaryOperation(Data const &arg_
2808        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2809        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2810        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2811          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
2812        }        }
2813        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2814          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2815        }        }
2816        // Compute a result for each tag        // Compute a result for each tag
2817        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2818        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2819          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2820          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2821          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);  
2822          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2823        }        }
2824    
# Line 2416  C_TensorBinaryOperation(Data const &arg_ Line 2834  C_TensorBinaryOperation(Data const &arg_
2834        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2835        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2836        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2837          res.requireWrite();
2838        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2839        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2840          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
2841          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2842          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2843            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2844            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2845            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2846            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2847            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2848          }          }
2849        }        }
2850    
2851      }      }
2852      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2853        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2854        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2855        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2856        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2857    
2858        int sampleNo_0,dataPointNo_0;        int sampleNo_0;
2859        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2860        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2861        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2862        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        const double *ptr_src = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2863          double ptr_1 = ptr_src[0];
2864          int size = size0 * numDataPointsPerSample_0;
2865          res.requireWrite();
2866          #pragma omp parallel for private(sampleNo_0) schedule(static)
2867        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2868          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {  //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2869            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0);
2870            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,0);
2871            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2872            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);  //          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2873            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2874            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
2875          }  //        }
2876        }        }
2877    
2878    
# Line 2466  C_TensorBinaryOperation(Data const &arg_ Line 2888  C_TensorBinaryOperation(Data const &arg_
2888        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2889        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2890        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2891          res.requireWrite();
2892        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2893        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2894          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2895          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2896          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2897            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2898            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2899            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2900            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2901            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2902          }          }
2903        }        }
# Line 2491  C_TensorBinaryOperation(Data const &arg_ Line 2914  C_TensorBinaryOperation(Data const &arg_
2914        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2915        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2916        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2917          res.requireWrite();
2918        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2919        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2920          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2921            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2922            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2923            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2924            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2925            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2926            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2927            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2928          }          }
2929        }        }
# Line 2521  Data Line 2945  Data
2945  C_TensorUnaryOperation(Data const &arg_0,  C_TensorUnaryOperation(Data const &arg_0,
2946                         UnaryFunction operation)                         UnaryFunction operation)
2947  {  {
2948      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2949      {
2950         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2951      }
2952      if (arg_0.isLazy())
2953      {
2954         throw DataException("Error - Operations not permitted on lazy data.");
2955      }
2956    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2957    Data arg_0_Z = Data(arg_0);    Data arg_0_Z = Data(arg_0);
2958    
2959    // Get rank and shape of inputs    // Get rank and shape of inputs
2960    int rank0 = arg_0_Z.getDataPointRank();    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
   DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();  
2961    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2962    
2963    // Declare output Data object    // Declare output Data object
# Line 2534  C_TensorUnaryOperation(Data const &arg_0 Line 2965  C_TensorUnaryOperation(Data const &arg_0
2965    
2966    if (arg_0_Z.isConstant()) {    if (arg_0_Z.isConstant()) {
2967      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2968      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2969      double *ptr_2 = &((res.getPointDataView().getData())[0]);      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2970      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2971    }    }
2972    else if (arg_0_Z.isTagged()) {    else if (arg_0_Z.isTagged()) {
# Line 2548  C_TensorUnaryOperation(Data const &arg_0 Line 2979  C_TensorUnaryOperation(Data const &arg_0
2979      res.tag();      res.tag();
2980      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2981    
     // Get the views  
     DataArrayView view_0 = tmp_0->getDefaultValue();  
     DataArrayView view_2 = tmp_2->getDefaultValue();  
2982      // Get the pointers to the actual data      // Get the pointers to the actual data
2983      double *ptr_0 = &((view_0.getData())[0]);      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2984      double *ptr_2 = &((view_2.getData())[0]);      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2985      // Compute a result for the default      // Compute a result for the default
2986      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2987      // Compute a result for each tag      // Compute a result for each tag
2988      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2989      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
2990      for (i=lookup_0.begin();i!=lookup_0.end();i++) {      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2991        tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());        tmp_2->addTag(i->first);
2992        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2993        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);  
2994        tensor_unary_operation(size0, ptr_0, ptr_2, operation);        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2995      }      }
2996    
# Line 2580  C_TensorUnaryOperation(Data const &arg_0 Line 3006  C_TensorUnaryOperation(Data const &arg_0
3006      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3007      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3008      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3009        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {      dataPointNo_0=0;
3010    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3011          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3012          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3013          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3014          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3015          tensor_unary_operation(size0, ptr_0, ptr_2, operation);          tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3016        }  //      }
3017      }      }
   
3018    }    }
3019    else {    else {
3020      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.2881

  ViewVC Help
Powered by ViewVC 1.1.26