/[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 2771 by jfenwick, Wed Nov 25 01:38:49 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /** \file Data.h */  /** \file Data.h */
16    
# Line 19  Line 18 
18  #define DATA_H  #define DATA_H
19  #include "system_dep.h"  #include "system_dep.h"
20    
21    #include "DataTypes.h"
22  #include "DataAbstract.h"  #include "DataAbstract.h"
23  #include "DataAlgorithm.h"  #include "DataAlgorithm.h"
24  #include "FunctionSpace.h"  #include "FunctionSpace.h"
# Line 26  Line 26 
26  #include "UnaryOp.h"  #include "UnaryOp.h"
27  #include "DataException.h"  #include "DataException.h"
28    
29    
30  extern "C" {  extern "C" {
31  #include "DataC.h"  #include "DataC.h"
32  /* #include "paso/Paso.h" doesn't belong in this file...causes trouble for BruceFactory.cpp */  //#include <omp.h>
33  }  }
34    
35    #ifdef _OPENMP
36    #include <omp.h>
37    #endif
38    
39  #include "esysmpi.h"  #include "esysmpi.h"
40  #include <string>  #include <string>
41  #include <algorithm>  #include <algorithm>
42    #include <sstream>
43    
44  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
45  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
46  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
 #include <boost/python/numeric.hpp>  
47    
48  namespace escript {  namespace escript {
49    
# Line 47  namespace escript { Line 52  namespace escript {
52  class DataConstant;  class DataConstant;
53  class DataTagged;  class DataTagged;
54  class DataExpanded;  class DataExpanded;
55    class DataLazy;
56    
57  /**  /**
58     \brief     \brief
59     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
60    
61     Description:     Description:
62     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
63     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
64     the object created for the lifetime of the object.     of the Data object.
65     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
66     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
67       Doing so will lead to invalid memory access.
68       This should not affect any methods exposed via boost::python.
69  */  */
70  class Data {  class Data {
71    
# Line 101  class Data { Line 108  class Data {
108         const FunctionSpace& what);         const FunctionSpace& what);
109    
110    /**    /**
111       \brief      \brief Copy Data from an existing vector
112       Constructor which copies data from a DataArrayView.    */
113    
      \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.  
   */  
114    ESCRIPT_DLL_API    ESCRIPT_DLL_API
115    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
116         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
117         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
118                     bool expanded=false);
119    
120    /**    /**
121       \brief       \brief
122       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
123    
124       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
125       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 128  class Data { Line 130  class Data {
130    */    */
131    ESCRIPT_DLL_API    ESCRIPT_DLL_API
132    Data(double value,    Data(double value,
133         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
134         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
135         bool expanded=false);         bool expanded=false);
136    
# Line 141  class Data { Line 143  class Data {
143    */    */
144    ESCRIPT_DLL_API    ESCRIPT_DLL_API
145    Data(const Data& inData,    Data(const Data& inData,
146         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
147    
148    /**    /**
149       \brief       \brief
150       Constructor which will create Tagged data if expanded is false.       Constructor which copies data from any object that can be treated like a python array/sequence.
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
     ==>*  
   */  
   ESCRIPT_DLL_API  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   ESCRIPT_DLL_API  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from any object that can be converted into  
      a python numarray.  
151    
152       \param value - Input - Input data.       \param value - Input - Input data.
153       \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 163  class Data {
163    /**    /**
164       \brief       \brief
165       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
166       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
167       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
168    
169       \param value - Input - Input data.       \param value - Input - Input data.
170       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
# Line 217  class Data { Line 182  class Data {
182         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
183         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
184         bool expanded=false);         bool expanded=false);
185    
186    
187    
188      /**
189        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
190        Once you have passed the pointer, do not delete it.
191      */
192      ESCRIPT_DLL_API
193      explicit Data(DataAbstract* underlyingdata);
194    
195      /**
196        \brief Create a Data based on the supplied DataAbstract
197      */
198      ESCRIPT_DLL_API
199      explicit Data(DataAbstract_ptr underlyingdata);
200    
201    /**    /**
202       \brief       \brief
203       Destructor       Destructor
# Line 225  class Data { Line 206  class Data {
206    ~Data();    ~Data();
207    
208    /**    /**
209       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
210    */    */
211    ESCRIPT_DLL_API    ESCRIPT_DLL_API
212    void    void
213    copy(const Data& other);    copy(const Data& other);
214    
215    /**    /**
216         \brief Return a pointer to a deep copy of this object.
217      */
218      ESCRIPT_DLL_API
219      Data
220      copySelf();
221    
222    
223      /**
224         \brief produce a delayed evaluation version of this Data.
225      */
226      ESCRIPT_DLL_API
227      Data
228      delay();
229    
230      /**
231         \brief convert the current data into lazy data.
232      */
233      ESCRIPT_DLL_API
234      void
235      delaySelf();
236    
237    
238      /**
239       Member access methods.       Member access methods.
240    */    */
241    
# Line 247  class Data { Line 250  class Data {
250    
251    /**    /**
252       \brief       \brief
253       Returns trueif the data object is protected against update       Returns true, if the data object is protected against update
254    
255    */    */
256    ESCRIPT_DLL_API    ESCRIPT_DLL_API
257    bool    bool
258    isProtected() const;    isProtected() const;
259    
260    /**  
261       \brief  /**
262       Return the values of a data point on this process     \brief
263    */     Return the value of a data point as a python tuple.
264    */
265    ESCRIPT_DLL_API    ESCRIPT_DLL_API
266    const boost::python::numeric::array    const boost::python::object
267    getValueOfDataPoint(int dataPointNo);    getValueOfDataPointAsTuple(int dataPointNo);
268    
269    /**    /**
270       \brief       \brief
# Line 272  class Data { Line 276  class Data {
276    
277    /**    /**
278       \brief       \brief
279       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
280    */    */
281    ESCRIPT_DLL_API    ESCRIPT_DLL_API
282    void    void
283    setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array&);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
284    
285    /**    /**
286       \brief       \brief
# Line 287  class Data { Line 291  class Data {
291    setValueOfDataPoint(int dataPointNo, const double);    setValueOfDataPoint(int dataPointNo, const double);
292    
293    /**    /**
294       \brief       \brief Return a data point across all processors as a python tuple.
      Return the value of the specified data-point across all processors  
295    */    */
296    ESCRIPT_DLL_API    ESCRIPT_DLL_API
297    const boost::python::numeric::array    const boost::python::object
298    getValueOfGlobalDataPoint(int procNo, int dataPointNo);    getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
299    
300    /**    /**
301       \brief       \brief
302       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
303    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
304    */    */
305    ESCRIPT_DLL_API    ESCRIPT_DLL_API
306    int    int
# Line 313  class Data { Line 314  class Data {
314    escriptDataC    escriptDataC
315    getDataC();    getDataC();
316    
317    
318    
319    /**    /**
320       \brief       \brief
321       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 324  class Data {
324    escriptDataC    escriptDataC
325    getDataC() const;    getDataC() const;
326    
   /**  
      \brief  
      Write the data as a string.  
   */  
   ESCRIPT_DLL_API  
   inline  
   std::string  
   toString() const  
   {  
     return m_data->toString();  
   }  
327    
328    /**    /**
329       \brief       \brief
330       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.  
331    */    */
332    ESCRIPT_DLL_API    ESCRIPT_DLL_API
333    inline    std::string
334    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
335    
336    /**    /**
337       \brief       \brief
# Line 360  class Data { Line 346  class Data {
346       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
347       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
348       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
349    */    */
350    ESCRIPT_DLL_API    ESCRIPT_DLL_API
351    void    void
352    tag();    tag();
353    
354    /**    /**
355        \brief If this data is lazy, then convert it to ready data.
356        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
357      */
358      ESCRIPT_DLL_API
359      void
360      resolve();
361    
362    
363      /**
364       \brief Ensures data is ready for write access.
365      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
366      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
367      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
368      Doing so might introduce additional sharing.
369      */
370      ESCRIPT_DLL_API
371      void
372      requireWrite();
373    
374      /**
375       \brief       \brief
376       Return true if this Data is expanded.       Return true if this Data is expanded.
377         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
378    */    */
379    ESCRIPT_DLL_API    ESCRIPT_DLL_API
380    bool    bool
# Line 376  class Data { Line 382  class Data {
382    
383    /**    /**
384       \brief       \brief
385         Return true if this Data is expanded or resolves to expanded.
386         That is, if it has a separate value for each datapoint in the sample.
387      */
388      ESCRIPT_DLL_API
389      bool
390      actsExpanded() const;
391      
392    
393      /**
394         \brief
395       Return true if this Data is tagged.       Return true if this Data is tagged.
396    */    */
397    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 391  class Data { Line 407  class Data {
407    isConstant() const;    isConstant() const;
408    
409    /**    /**
410         \brief Return true if this Data is lazy.
411      */
412      ESCRIPT_DLL_API
413      bool
414      isLazy() const;
415    
416      /**
417         \brief Return true if this data is ready.
418      */
419      ESCRIPT_DLL_API
420      bool
421      isReady() const;
422    
423      /**
424       \brief       \brief
425       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
426    contains datapoints.
427    */    */
428    ESCRIPT_DLL_API    ESCRIPT_DLL_API
429    bool    bool
# Line 424  class Data { Line 455  class Data {
455    */    */
456    ESCRIPT_DLL_API    ESCRIPT_DLL_API
457    inline    inline
458    const AbstractDomain&  //   const AbstractDomain&
459      const_Domain_ptr
460    getDomain() const    getDomain() const
461    {    {
462       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
463    }    }
464    
465    
466      /**
467         \brief
468         Return the domain.
469         TODO: For internal use only.   This should be removed.
470      */
471      ESCRIPT_DLL_API
472      inline
473    //   const AbstractDomain&
474      Domain_ptr
475      getDomainPython() const
476      {
477         return getFunctionSpace().getDomainPython();
478      }
479    
480    /**    /**
481       \brief       \brief
482       Return a copy of the domain.       Return a copy of the domain.
# Line 444  class Data { Line 491  class Data {
491    */    */
492    ESCRIPT_DLL_API    ESCRIPT_DLL_API
493    inline    inline
494    int    unsigned int
495    getDataPointRank() const    getDataPointRank() const
496    {    {
497      return m_data->getPointDataView().getRank();      return m_data->getRank();
498    }    }
499    
500    /**    /**
# Line 484  class Data { Line 531  class Data {
531    {    {
532      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
533    }    }
534    
535    
536      /**
537        \brief
538        Return the number of values in the shape for this object.
539      */
540      ESCRIPT_DLL_API
541      int
542      getNoValues() const
543      {
544        return m_data->getNoValues();
545      }
546    
547    
548    /**    /**
549       \brief       \brief
550       dumps the object into a netCDF file       dumps the object into a netCDF file
# Line 491  class Data { Line 552  class Data {
552    ESCRIPT_DLL_API    ESCRIPT_DLL_API
553    void    void
554    dump(const std::string fileName) const;    dump(const std::string fileName) const;
555    
556     /**
557      \brief returns the values of the object as a list of tuples (one for each datapoint).
558    
559      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
560    If false, the result is a list of scalars [1, 2, ...]
561     */
562      ESCRIPT_DLL_API
563      const boost::python::object
564      toListOfTuples(bool scalarastuple=true);
565    
566    
567     /**
568        \brief
569        Return the sample data for the given sample no. This is not the
570        preferred interface but is provided for use by C code.
571        The bufferg parameter is only required for LazyData.
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.
611       NOTE: Construction of the DataArrayView is a relatively expensive       \param sampleNo - Input -
612       operation.       \param dataPointNo - Input -
613      */
614      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 -       \param sampleNo - Input -
622       \param dataPointNo - Input -       \param dataPointNo - Input -
623    */    */
624    ESCRIPT_DLL_API    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
812      boost::python::object
813      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    ESCRIPT_DLL_API
821    boost::python::numeric::array    boost::python::object
822    integrate() const;    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    
1760    };
1761    
1762    }   // end namespace escript
1763    
1764    
1765    // No, this is not supposed to be at the top of the file
1766    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1767    // so that I can dynamic cast between them below.
1768    #include "DataReady.h"
1769    #include "DataLazy.h"
1770    
1771    namespace escript
1772    {
1773    
1774    inline
1775    const DataReady*
1776    Data::getReady() const
1777    {
1778       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1779       EsysAssert((dr!=0), "Error - casting to DataReady.");
1780       return dr;
1781    }
1782    
1783    inline
1784    DataReady*
1785    Data::getReady()
1786    {
1787       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1788       EsysAssert((dr!=0), "Error - casting to DataReady.");
1789       return dr;
1790    }
1791    
1792    // Be wary of using this for local operations since it (temporarily) increases reference count.
1793    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1794    // getReady() instead
1795    inline
1796    DataReady_ptr
1797    Data::getReadyPtr()
1798    {
1799       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1800       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1801       return dr;
1802    }
1803    
1804    
1805    inline
1806    const_DataReady_ptr
1807    Data::getReadyPtr() const
1808    {
1809       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1810       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1811       return dr;
1812  }  }
1813    
1814    inline
1815    DataAbstract::ValueType::value_type*
1816    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1817    {
1818       if (isLazy())
1819       {
1820        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1821       }
1822       return getReady()->getSampleDataRW(sampleNo);
1823    }
1824    
1825    inline
1826    const DataAbstract::ValueType::value_type*
1827    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo)
1828    {
1829       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1830       if (l!=0)
1831       {
1832        size_t offset=0;
1833        const DataTypes::ValueType* res=l->resolveSample(sampleNo,offset);
1834        return &((*res)[offset]);
1835       }
1836       return getReady()->getSampleDataRO(sampleNo);
1837    }
1838    
1839    
1840    
1841    /**
1842       Modify a filename for MPI parallel output to multiple files
1843    */
1844    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1845    
1846  /**  /**
1847     Binary Data object operators.     Binary Data object operators.
1848  */  */
# Line 1519  ESCRIPT_DLL_API std::ostream& operator<< Line 1954  ESCRIPT_DLL_API std::ostream& operator<<
1954  /**  /**
1955    \brief    \brief
1956    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1957    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1958    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1959    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1960    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1961  */  */
1962  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1963  Data  Data
1964  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1965                       Data& arg1,                       Data& arg_1,
1966                       int axis_offset=0,                       int axis_offset=0,
1967                       int transpose=0);                       int transpose=0);
1968    
   
 /**  
   \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);  
   
1969  /**  /**
1970    \brief    \brief
1971    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 1979  Data::binaryOp(const Data& right,
1979  {  {
1980     //     //
1981     // 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
1982     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1983       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.");
1984     }     }
1985    
1986       if (isLazy() || right.isLazy())
1987       {
1988         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1989       }
1990     //     //
1991     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1992     Data tempRight(right);     Data tempRight(right);
1993    
1994     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1995       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1996         //         //
# Line 1582  Data::binaryOp(const Data& right, Line 2000  Data::binaryOp(const Data& right,
2000         //         //
2001         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2002         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2003         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2004           set_m_data(tempLeft.m_data);
2005       }       }
2006     }     }
2007     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1598  Data::binaryOp(const Data& right, Line 2017  Data::binaryOp(const Data& right,
2017       // of any data type       // of any data type
2018       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2019       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2020       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2021     } else if (isTagged()) {     } else if (isTagged()) {
2022       //       //
2023       // 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 2065  Data::algorithm(BinaryFunction operation
2065      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2066      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2067      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2068      } else if (isEmpty()) {
2069        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2070      } else if (isLazy()) {
2071        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2072      } else {
2073        throw DataException("Error - Data encapsulates an unknown type.");
2074    }    }
   return 0;  
2075  }  }
2076    
2077  /**  /**
# Line 1663  inline Line 2087  inline
2087  Data  Data
2088  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2089  {  {
2090    if (isExpanded()) {    if (isEmpty()) {
2091      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2092      }
2093      else if (isExpanded()) {
2094        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2095      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2096      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2097      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2098      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2099      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2100      return result;      return result;
2101    } else if (isTagged()) {    }
2102      else if (isTagged()) {
2103      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());  
2104      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2105      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2106        defval[0]=0;
2107        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2108      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2109      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2110    } else if (isConstant()) {    }
2111      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2112        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2113      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2114      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2115      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2116      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2117      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2118      return result;      return result;
2119      } else if (isLazy()) {
2120        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2121      } else {
2122        throw DataException("Error - Data encapsulates an unknown type.");
2123    }    }
   Data falseRetVal; // to keep compiler quiet  
   return falseRetVal;  
2124  }  }
2125    
2126    /**
2127      \brief
2128      Compute a tensor operation with two Data objects
2129      \param arg_0 - Input - Data object
2130      \param arg_1 - Input - Data object
2131      \param operation - Input - Binary op functor
2132    */
2133  template <typename BinaryFunction>  template <typename BinaryFunction>
2134    inline
2135  Data  Data
2136  C_TensorBinaryOperation(Data const &arg_0,  C_TensorBinaryOperation(Data const &arg_0,
2137                          Data const &arg_1,                          Data const &arg_1,
2138                          BinaryFunction operation)                          BinaryFunction operation)
2139  {  {
2140      if (arg_0.isEmpty() || arg_1.isEmpty())
2141      {
2142         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2143      }
2144      if (arg_0.isLazy() || arg_1.isLazy())
2145      {
2146         throw DataException("Error - Operations not permitted on lazy data.");
2147      }
2148    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2149    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
2150    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
# Line 1730  C_TensorBinaryOperation(Data const &arg_ Line 2166  C_TensorBinaryOperation(Data const &arg_
2166    // Get rank and shape of inputs    // Get rank and shape of inputs
2167    int rank0 = arg_0_Z.getDataPointRank();    int rank0 = arg_0_Z.getDataPointRank();
2168    int rank1 = arg_1_Z.getDataPointRank();    int rank1 = arg_1_Z.getDataPointRank();
2169    DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();    DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2170    DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();    DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2171    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2172    int size1 = arg_1_Z.getDataPointSize();    int size1 = arg_1_Z.getDataPointSize();
   
2173    // Declare output Data object    // Declare output Data object
2174    Data res;    Data res;
2175    
2176    if (shape0 == shape1) {    if (shape0 == shape1) {
   
2177      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2178        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2179        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2180        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2181        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2182    
2183        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2184      }      }
2185      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 2197  C_TensorBinaryOperation(Data const &arg_
2197    
2198        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2199        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2200        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2201        // Get the views  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2202        // Get the pointers to the actual data        // Get the pointers to the actual data
2203        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2204        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2205    
2206        // Compute a result for the default        // Compute a result for the default
2207        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2208        // Compute a result for each tag        // Compute a result for each tag
2209        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2210        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
2211        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2212          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2213          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2214          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2215          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2216          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2217        }        }
2218    
2219      }      }
2220      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
   
2221        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2222        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2223        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 2227  C_TensorBinaryOperation(Data const &arg_
2227        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2228        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2229        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2230          res.requireWrite();
2231        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2232        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2233          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2234            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2235            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2236            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2237            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2238            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2239            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2240          }          }
2241        }        }
2242    
2243      }      }
2244      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
   
2245        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2246        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2247    
# Line 1823  C_TensorBinaryOperation(Data const &arg_ Line 2255  C_TensorBinaryOperation(Data const &arg_
2255    
2256        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2257        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2258        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);  
2259        // 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();  
2260        // Get the pointers to the actual data        // Get the pointers to the actual data
2261        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2262        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2263        // Compute a result for the default        // Compute a result for the default
2264        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2265        // Compute a result for each tag        // Compute a result for each tag
2266        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2267        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
2268        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2269          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2270          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2271          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);  
2272          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2273        }        }
2274    
2275      }      }
2276      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
   
2277        // Borrow DataTagged input from Data object        // Borrow DataTagged input from Data object
2278        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2279    
# Line 1858  C_TensorBinaryOperation(Data const &arg_ Line 2285  C_TensorBinaryOperation(Data const &arg_
2285        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2286        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2287    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2288        // Get the pointers to the actual data        // Get the pointers to the actual data
2289        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2290        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2291        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2292    
2293        // Compute a result for the default        // Compute a result for the default
2294        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2295        // Merge the tags        // Merge the tags
# Line 1873  C_TensorBinaryOperation(Data const &arg_ Line 2297  C_TensorBinaryOperation(Data const &arg_
2297        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2298        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2299        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2300          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape          tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2301        }        }
2302        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2303          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2304        }        }
2305        // Compute a result for each tag        // Compute a result for each tag
2306        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2307        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2308          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);  
2309          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2310          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2311          double *ptr_0 = &view_0.getData(0);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2312          double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2313          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2314        }        }
2315    
2316      }      }
2317      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
   
2318        // 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
2319        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2320        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 2324  C_TensorBinaryOperation(Data const &arg_
2324        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2325        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2326        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2327          res.requireWrite();
2328        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2329        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2330          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
2331          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2332          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2333            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2334            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2335            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2336            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2337            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2338          }          }
2339        }        }
2340    
2341      }      }
2342      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2343        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2344        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2345        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 2349  C_TensorBinaryOperation(Data const &arg_
2349        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2350        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2351        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2352          res.requireWrite();
2353        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2354        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2355          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2356            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2357            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2358            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);  
2359            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2360            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2361              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2362    
2363    
2364            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2365          }          }
2366        }        }
2367    
2368      }      }
2369      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
   
2370        // 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
2371        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2372        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 2376  C_TensorBinaryOperation(Data const &arg_
2376        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2377        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2378        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2379          res.requireWrite();
2380        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2381        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2382          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2383          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2384          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2385            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2386            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2387            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2388            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2389            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2390          }          }
2391        }        }
2392    
2393      }      }
2394      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
   
2395        // 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
2396        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2397        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 2401  C_TensorBinaryOperation(Data const &arg_
2401        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2402        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2403        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2404          res.requireWrite();
2405        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2406        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2407          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        dataPointNo_0=0;
2408    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2409            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2410            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2411            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2412            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2413            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2414            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2415            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);            tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2416          }  //       }
2417        }        }
2418    
2419      }      }
# Line 1995  C_TensorBinaryOperation(Data const &arg_ Line 2422  C_TensorBinaryOperation(Data const &arg_
2422      }      }
2423    
2424    } else if (0 == rank0) {    } else if (0 == rank0) {
   
2425      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2426        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2427        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2428        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2429        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2430        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2431      }      }
2432      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 2444  C_TensorBinaryOperation(Data const &arg_
2444    
2445        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2446        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2447        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2448        // Get the views  
2449        DataArrayView view_1 = tmp_1->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2450        DataArrayView view_2 = tmp_2->getDefaultValue();        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2451        // Get the pointers to the actual data  
       double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2452        // Compute a result for the default        // Compute a result for the default
2453        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2454        // Compute a result for each tag        // Compute a result for each tag
2455        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2456        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
2457        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2458          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2459          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2460          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);  
2461          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2462        }        }
2463    
# Line 2051  C_TensorBinaryOperation(Data const &arg_ Line 2473  C_TensorBinaryOperation(Data const &arg_
2473        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2474        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2475        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2476          res.requireWrite();
2477        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2478        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2479          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2480            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2481            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2482            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2483            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2484            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2485            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2486    
2487          }          }
# Line 2080  C_TensorBinaryOperation(Data const &arg_ Line 2503  C_TensorBinaryOperation(Data const &arg_
2503    
2504        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2505        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2506        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2507        // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2508        // Get the pointers to the actual data        // Get the pointers to the actual data
2509        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2510        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2511    
2512    
2513        // Compute a result for the default        // Compute a result for the default
2514        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2515        // Compute a result for each tag        // Compute a result for each tag
2516        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2517        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
2518        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2519          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2520          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2521          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2522          double *ptr_0 = &view_0.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2523          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2524        }        }
2525    
# Line 2115  C_TensorBinaryOperation(Data const &arg_ Line 2537  C_TensorBinaryOperation(Data const &arg_
2537        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2538        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2539    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2540        // Get the pointers to the actual data        // Get the pointers to the actual data
2541        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2542        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2543        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2544    
2545        // Compute a result for the default        // Compute a result for the default
2546        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2547        // Merge the tags        // Merge the tags
# Line 2130  C_TensorBinaryOperation(Data const &arg_ Line 2549  C_TensorBinaryOperation(Data const &arg_
2549        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2550        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2551        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2552          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
2553        }        }
2554        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2555          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2556        }        }
2557        // Compute a result for each tag        // Compute a result for each tag
2558        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2559        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2560          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2561          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2562          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2563          double *ptr_0 = &view_0.getData(0);  
         double *ptr_1 = &view_1.getData(0);  
         double *ptr_2 = &view_2.getData(0);  
2564          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2565        }        }
2566    
# Line 2159  C_TensorBinaryOperation(Data const &arg_ Line 2576  C_TensorBinaryOperation(Data const &arg_
2576        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2577        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2578        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2579          res.requireWrite();
2580        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2581        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2582          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
2583          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2584          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2585            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2586            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2587            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2588            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2589            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2590          }          }
2591        }        }
2592    
2593      }      }
2594      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2595        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2596        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2597        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 2601  C_TensorBinaryOperation(Data const &arg_
2601        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2602        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2603        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2604          res.requireWrite();
2605        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2606        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2607          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2608            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2609            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2610            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2611            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2612            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2613            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2614          }          }
2615        }        }
# Line 2209  C_TensorBinaryOperation(Data const &arg_ Line 2627  C_TensorBinaryOperation(Data const &arg_
2627        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2628        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2629        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2630          res.requireWrite();
2631        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2632        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2633          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2634          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2635          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2636            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2637            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2638            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2639            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2640            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2641          }          }
2642        }        }
# Line 2234  C_TensorBinaryOperation(Data const &arg_ Line 2653  C_TensorBinaryOperation(Data const &arg_
2653        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2654        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2655        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2656          res.requireWrite();
2657        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2658        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2659          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2660            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2661            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2662            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2663            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2664            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2665            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2666            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2667          }          }
2668        }        }
# Line 2253  C_TensorBinaryOperation(Data const &arg_ Line 2673  C_TensorBinaryOperation(Data const &arg_
2673      }      }
2674    
2675    } else if (0 == rank1) {    } else if (0 == rank1) {
   
2676      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {      if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2677        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2678        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2679        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2680        double *ptr_2 = &((res.getPointDataView().getData())[0]);        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2681        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2682      }      }
2683      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 2695  C_TensorBinaryOperation(Data const &arg_
2695    
2696        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2697        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2698        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2699        // Get the views  
2700        DataArrayView view_1 = tmp_1->getDefaultValue();        //Get the pointers to the actual data
2701        DataArrayView view_2 = tmp_2->getDefaultValue();        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2702        // Get the pointers to the actual data        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2703        double *ptr_1 = &((view_1.getData())[0]);  
       double *ptr_2 = &((view_2.getData())[0]);  
2704        // Compute a result for the default        // Compute a result for the default
2705        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2706        // Compute a result for each tag        // Compute a result for each tag
2707        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2708        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
2709        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2710          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2711          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2712          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);  
2713          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2714        }        }
   
2715      }      }
2716      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2717    
# Line 2309  C_TensorBinaryOperation(Data const &arg_ Line 2724  C_TensorBinaryOperation(Data const &arg_
2724        int numSamples_1 = arg_1_Z.getNumSamples();        int numSamples_1 = arg_1_Z.getNumSamples();
2725        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2726        int offset_0 = tmp_0->getPointOffset(0,0);        int offset_0 = tmp_0->getPointOffset(0,0);
2727          res.requireWrite();
2728        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2729        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2730          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2731            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2732            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2733            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2734            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2735            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2736            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2737          }          }
2738        }        }
# Line 2337  C_TensorBinaryOperation(Data const &arg_ Line 2753  C_TensorBinaryOperation(Data const &arg_
2753    
2754        // Prepare offset into DataConstant        // Prepare offset into DataConstant
2755        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2756        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();  
2757        // Get the pointers to the actual data        // Get the pointers to the actual data
2758        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2759        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2760        // Compute a result for the default        // Compute a result for the default
2761        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2762        // Compute a result for each tag        // Compute a result for each tag
2763        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2764        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
2765        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2766          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2767          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2768          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);  
2769          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2770        }        }
2771    
# Line 2372  C_TensorBinaryOperation(Data const &arg_ Line 2783  C_TensorBinaryOperation(Data const &arg_
2783        res.tag();        // DataTagged output        res.tag();        // DataTagged output
2784        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2785    
       // Get the views  
       DataArrayView view_0 = tmp_0->getDefaultValue();  
       DataArrayView view_1 = tmp_1->getDefaultValue();  
       DataArrayView view_2 = tmp_2->getDefaultValue();  
2786        // Get the pointers to the actual data        // Get the pointers to the actual data
2787        double *ptr_0 = &((view_0.getData())[0]);        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2788        double *ptr_1 = &((view_1.getData())[0]);        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2789        double *ptr_2 = &((view_2.getData())[0]);        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2790    
2791        // Compute a result for the default        // Compute a result for the default
2792        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2793        // Merge the tags        // Merge the tags
# Line 2387  C_TensorBinaryOperation(Data const &arg_ Line 2795  C_TensorBinaryOperation(Data const &arg_
2795        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2796        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2797        for (i=lookup_0.begin();i!=lookup_0.end();i++) {        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2798          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
2799        }        }
2800        for (i=lookup_1.begin();i!=lookup_1.end();i++) {        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2801          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());          tmp_2->addTag(i->first);
2802        }        }
2803        // Compute a result for each tag        // Compute a result for each tag
2804        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2805        for (i=lookup_2.begin();i!=lookup_2.end();i++) {        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2806          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2807          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2808          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);  
2809          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2810        }        }
2811    
# Line 2416  C_TensorBinaryOperation(Data const &arg_ Line 2821  C_TensorBinaryOperation(Data const &arg_
2821        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2822        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2823        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2824          res.requireWrite();
2825        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2826        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2827          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
2828          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2829          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2830            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2831            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2832            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2833            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2834            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2835          }          }
2836        }        }
2837    
2838      }      }
2839      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
   
2840        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output        res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2841        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2842        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
# Line 2441  C_TensorBinaryOperation(Data const &arg_ Line 2846  C_TensorBinaryOperation(Data const &arg_
2846        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2847        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2848        int offset_1 = tmp_1->getPointOffset(0,0);        int offset_1 = tmp_1->getPointOffset(0,0);
2849          res.requireWrite();
2850        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2851        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2852          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2853            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2854            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2855            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2856            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2857            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2858            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2859          }          }
2860        }        }
# Line 2466  C_TensorBinaryOperation(Data const &arg_ Line 2872  C_TensorBinaryOperation(Data const &arg_
2872        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2873        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2874        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2875          res.requireWrite();
2876        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2877        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2878          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2879          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2880          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2881            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2882            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2883            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2884            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2885            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2886          }          }
2887        }        }
# Line 2491  C_TensorBinaryOperation(Data const &arg_ Line 2898  C_TensorBinaryOperation(Data const &arg_
2898        int sampleNo_0,dataPointNo_0;        int sampleNo_0,dataPointNo_0;
2899        int numSamples_0 = arg_0_Z.getNumSamples();        int numSamples_0 = arg_0_Z.getNumSamples();
2900        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2901          res.requireWrite();
2902        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2903        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2904          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2905            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2906            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2907            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2908            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2909            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2910            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2911            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2912          }          }
2913        }        }
# Line 2521  Data Line 2929  Data
2929  C_TensorUnaryOperation(Data const &arg_0,  C_TensorUnaryOperation(Data const &arg_0,
2930                         UnaryFunction operation)                         UnaryFunction operation)
2931  {  {
2932      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2933      {
2934         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2935      }
2936      if (arg_0.isLazy())
2937      {
2938         throw DataException("Error - Operations not permitted on lazy data.");
2939      }
2940    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2941    Data arg_0_Z = Data(arg_0);    Data arg_0_Z = Data(arg_0);
2942    
2943    // Get rank and shape of inputs    // Get rank and shape of inputs
2944    int rank0 = arg_0_Z.getDataPointRank();    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
   DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();  
2945    int size0 = arg_0_Z.getDataPointSize();    int size0 = arg_0_Z.getDataPointSize();
2946    
2947    // Declare output Data object    // Declare output Data object
# Line 2534  C_TensorUnaryOperation(Data const &arg_0 Line 2949  C_TensorUnaryOperation(Data const &arg_0
2949    
2950    if (arg_0_Z.isConstant()) {    if (arg_0_Z.isConstant()) {
2951      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2952      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2953      double *ptr_2 = &((res.getPointDataView().getData())[0]);      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2954      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2955    }    }
2956    else if (arg_0_Z.isTagged()) {    else if (arg_0_Z.isTagged()) {
# Line 2548  C_TensorUnaryOperation(Data const &arg_0 Line 2963  C_TensorUnaryOperation(Data const &arg_0
2963      res.tag();      res.tag();
2964      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2965    
     // Get the views  
     DataArrayView view_0 = tmp_0->getDefaultValue();  
     DataArrayView view_2 = tmp_2->getDefaultValue();  
2966      // Get the pointers to the actual data      // Get the pointers to the actual data
2967      double *ptr_0 = &((view_0.getData())[0]);      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2968      double *ptr_2 = &((view_2.getData())[0]);      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2969      // Compute a result for the default      // Compute a result for the default
2970      tensor_unary_operation(size0, ptr_0, ptr_2, operation);      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2971      // Compute a result for each tag      // Compute a result for each tag
2972      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2973      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
2974      for (i=lookup_0.begin();i!=lookup_0.end();i++) {      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2975        tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());        tmp_2->addTag(i->first);
2976        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2977        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);  
2978        tensor_unary_operation(size0, ptr_0, ptr_2, operation);        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2979      }      }
2980    
# Line 2580  C_TensorUnaryOperation(Data const &arg_0 Line 2990  C_TensorUnaryOperation(Data const &arg_0
2990      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2991      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2992      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2993        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {      dataPointNo_0=0;
2994    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2995          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2996          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2997          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2998          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2999          tensor_unary_operation(size0, ptr_0, ptr_2, operation);          tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3000        }  //      }
3001      }      }
   
3002    }    }
3003    else {    else {
3004      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.2771

  ViewVC Help
Powered by ViewVC 1.1.26