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

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

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

revision 922 by gross, Fri Jan 5 04:23:05 2007 UTC revision 2770 by jfenwick, Wed Nov 25 01:24:51 2009 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2009 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
8   *  Licensed under the Open Software License version 3.0    *  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Licensed under the Open Software License version 3.0
10   *                                                          *  * http://www.opensource.org/licenses/osl-3.0.php
11   ************************************************************  *
12  */  *******************************************************/
13    
14    
15  /** \file Data.h */  /** \file Data.h */
16    
# Line 17  Line 18 
18  #define DATA_H  #define DATA_H
19  #include "system_dep.h"  #include "system_dep.h"
20    
21    #include "DataTypes.h"
22  #include "DataAbstract.h"  #include "DataAbstract.h"
23  #include "DataAlgorithm.h"  #include "DataAlgorithm.h"
24  #include "FunctionSpace.h"  #include "FunctionSpace.h"
# Line 24  Line 26 
26  #include "UnaryOp.h"  #include "UnaryOp.h"
27  #include "DataException.h"  #include "DataException.h"
28    
29    
30  extern "C" {  extern "C" {
31  #include "DataC.h"  #include "DataC.h"
32  #include "paso/Paso.h"  //#include <omp.h>
33  }  }
34    
35  #ifndef PASO_MPI  #include "esysmpi.h"
 #define MPI_Comm long  
 #endif  
   
36  #include <string>  #include <string>
37  #include <algorithm>  #include <algorithm>
38    #include <sstream>
39    
40  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
41  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
42  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
 #include <boost/python/numeric.hpp>  
43    
44  namespace escript {  namespace escript {
45    
# Line 48  namespace escript { Line 48  namespace escript {
48  class DataConstant;  class DataConstant;
49  class DataTagged;  class DataTagged;
50  class DataExpanded;  class DataExpanded;
51    class DataLazy;
52    
53  /**  /**
54     \brief     \brief
55     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
56    
57     Description:     Description:
58     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
59     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
60     the object created for the lifetime of the object.     of the Data object.
61     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
62     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
63       Doing so will lead to invalid memory access.
64       This should not affect any methods exposed via boost::python.
65  */  */
66  class Data {  class Data {
67    
# Line 102  class Data { Line 104  class Data {
104         const FunctionSpace& what);         const FunctionSpace& what);
105    
106    /**    /**
107       \brief      \brief Copy Data from an existing vector
108       Constructor which copies data from a DataArrayView.    */
109    
      \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.  
   */  
110    ESCRIPT_DLL_API    ESCRIPT_DLL_API
111    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
112         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
113         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
114                     bool expanded=false);
115    
116    /**    /**
117       \brief       \brief
118       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
119    
120       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
121       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 129  class Data { Line 126  class Data {
126    */    */
127    ESCRIPT_DLL_API    ESCRIPT_DLL_API
128    Data(double value,    Data(double value,
129         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
130         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
131         bool expanded=false);         bool expanded=false);
132    
# Line 142  class Data { Line 139  class Data {
139    */    */
140    ESCRIPT_DLL_API    ESCRIPT_DLL_API
141    Data(const Data& inData,    Data(const Data& inData,
142         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Constructor which will create Tagged data if expanded is false.  
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
     ==>*  
   */  
   ESCRIPT_DLL_API  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   ESCRIPT_DLL_API  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
143    
144    /**    /**
145       \brief       \brief
146       Constructor which copies data from any object that can be converted into       Constructor which copies data from any object that can be treated like a python array/sequence.
      a python numarray.  
147    
148       \param value - Input - Input data.       \param value - Input - Input data.
149       \param what - Input - A description of what this data represents.       \param what - Input - A description of what this data represents.
# Line 199  class Data { Line 159  class Data {
159    /**    /**
160       \brief       \brief
161       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
162       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
163       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
164    
165       \param value - Input - Input data.       \param value - Input - Input data.
166       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
# Line 214  class Data { Line 174  class Data {
174       Constructor which creates a DataConstant of "shape" with constant value.       Constructor which creates a DataConstant of "shape" with constant value.
175    */    */
176    ESCRIPT_DLL_API    ESCRIPT_DLL_API
177    Data(double value,    Data(double value,
178         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
179         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
180         bool expanded=false);         bool expanded=false);
181    
182    
183    
184      /**
185        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
186        Once you have passed the pointer, do not delete it.
187      */
188      ESCRIPT_DLL_API
189      explicit Data(DataAbstract* underlyingdata);
190    
191      /**
192        \brief Create a Data based on the supplied DataAbstract
193      */
194      ESCRIPT_DLL_API
195      explicit Data(DataAbstract_ptr underlyingdata);
196    
197    /**    /**
198       \brief       \brief
199       Destructor       Destructor
# Line 226  class Data { Line 202  class Data {
202    ~Data();    ~Data();
203    
204    /**    /**
205       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
206    */    */
207    ESCRIPT_DLL_API    ESCRIPT_DLL_API
208    void    void
209    copy(const Data& other);    copy(const Data& other);
210    
211    /**    /**
212         \brief Return a pointer to a deep copy of this object.
213      */
214      ESCRIPT_DLL_API
215      Data
216      copySelf();
217    
218    
219      /**
220         \brief produce a delayed evaluation version of this Data.
221      */
222      ESCRIPT_DLL_API
223      Data
224      delay();
225    
226      /**
227         \brief convert the current data into lazy data.
228      */
229      ESCRIPT_DLL_API
230      void
231      delaySelf();
232    
233    
234      /**
235       Member access methods.       Member access methods.
236    */    */
237    
238    /**    /**
239       \brief       \brief
240       switches on update protection       switches on update protection
241    
242    */    */
243    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 248  class Data { Line 246  class Data {
246    
247    /**    /**
248       \brief       \brief
249       Returns trueif the data object is protected against update       Returns true, if the data object is protected against update
250    
251    */    */
252    ESCRIPT_DLL_API    ESCRIPT_DLL_API
253    bool    bool
254    isProtected() const;    isProtected() const;
   /**  
      \brief  
      Return the values of all data-points as a single python numarray object.  
   */  
   ESCRIPT_DLL_API  
   const boost::python::numeric::array  
   convertToNumArray();  
255    
256    /**  
257       \brief  /**
258       Fills the expanded Data object from values of a python numarray object.     \brief
259    */     Return the value of a data point as a python tuple.
260    */
261    ESCRIPT_DLL_API    ESCRIPT_DLL_API
262    void    const boost::python::object
263    fillFromNumArray(const boost::python::numeric::array);    getValueOfDataPointAsTuple(int dataPointNo);
264    
265    /**    /**
266       \brief       \brief
267       Return the values of a data point on this process       sets the values of a data-point from a python object on this process
268    */    */
269    ESCRIPT_DLL_API    ESCRIPT_DLL_API
270    const boost::python::numeric::array    void
271    getValueOfDataPoint(int dataPointNo);    setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
272    
273    /**    /**
274       \brief       \brief
275       sets the values of a data-point on this process       sets the values of a data-point from a array-like object on this process
276    */    */
277    ESCRIPT_DLL_API    ESCRIPT_DLL_API
278    void    void
279    setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array);    setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
280    
281    /**    /**
282       \brief       \brief
# Line 295  class Data { Line 287  class Data {
287    setValueOfDataPoint(int dataPointNo, const double);    setValueOfDataPoint(int dataPointNo, const double);
288    
289    /**    /**
290       \brief       \brief Return a data point across all processors as a python tuple.
      Return the value of the specified data-point across all processors  
291    */    */
292    ESCRIPT_DLL_API    ESCRIPT_DLL_API
293    const boost::python::numeric::array    const boost::python::object
294    getValueOfGlobalDataPoint(int procNo, int dataPointNo);    getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
295    
296    /**    /**
297       \brief       \brief
298       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
299    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
300    */    */
301    ESCRIPT_DLL_API    ESCRIPT_DLL_API
302    int    int
# Line 321  class Data { Line 310  class Data {
310    escriptDataC    escriptDataC
311    getDataC();    getDataC();
312    
313    
314    
315    /**    /**
316       \brief       \brief
317       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 329  class Data { Line 320  class Data {
320    escriptDataC    escriptDataC
321    getDataC() const;    getDataC() const;
322    
   /**  
      \brief  
      Write the data as a string.  
   */  
   ESCRIPT_DLL_API  
   inline  
   std::string  
   toString() const  
   {  
     return m_data->toString();  
   }  
323    
324    /**    /**
325       \brief       \brief
326       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.  
327    */    */
328    ESCRIPT_DLL_API    ESCRIPT_DLL_API
329    inline    std::string
330    const DataArrayView&    toString() const;
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
331    
332    /**    /**
333       \brief       \brief
# Line 368  class Data { Line 342  class Data {
342       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
343       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
344       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
345    */    */
346    ESCRIPT_DLL_API    ESCRIPT_DLL_API
347    void    void
348    tag();    tag();
349    
350    /**    /**
351        \brief If this data is lazy, then convert it to ready data.
352        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
353      */
354      ESCRIPT_DLL_API
355      void
356      resolve();
357    
358    
359      /**
360       \brief Ensures data is ready for write access.
361      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
362      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
363      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
364      Doing so might introduce additional sharing.
365      */
366      ESCRIPT_DLL_API
367      void
368      requireWrite();
369    
370      /**
371       \brief       \brief
372       Return true if this Data is expanded.       Return true if this Data is expanded.
373         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
374    */    */
375    ESCRIPT_DLL_API    ESCRIPT_DLL_API
376    bool    bool
# Line 384  class Data { Line 378  class Data {
378    
379    /**    /**
380       \brief       \brief
381         Return true if this Data is expanded or resolves to expanded.
382         That is, if it has a separate value for each datapoint in the sample.
383      */
384      ESCRIPT_DLL_API
385      bool
386      actsExpanded() const;
387      
388    
389      /**
390         \brief
391       Return true if this Data is tagged.       Return true if this Data is tagged.
392    */    */
393    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 399  class Data { Line 403  class Data {
403    isConstant() const;    isConstant() const;
404    
405    /**    /**
406         \brief Return true if this Data is lazy.
407      */
408      ESCRIPT_DLL_API
409      bool
410      isLazy() const;
411    
412      /**
413         \brief Return true if this data is ready.
414      */
415      ESCRIPT_DLL_API
416      bool
417      isReady() const;
418    
419      /**
420       \brief       \brief
421       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
422    contains datapoints.
423    */    */
424    ESCRIPT_DLL_API    ESCRIPT_DLL_API
425    bool    bool
# Line 432  class Data { Line 451  class Data {
451    */    */
452    ESCRIPT_DLL_API    ESCRIPT_DLL_API
453    inline    inline
454    const AbstractDomain&  //   const AbstractDomain&
455      const_Domain_ptr
456    getDomain() const    getDomain() const
457    {    {
458       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
459    }    }
460    
461    
462      /**
463         \brief
464         Return the domain.
465         TODO: For internal use only.   This should be removed.
466      */
467      ESCRIPT_DLL_API
468      inline
469    //   const AbstractDomain&
470      Domain_ptr
471      getDomainPython() const
472      {
473         return getFunctionSpace().getDomainPython();
474      }
475    
476    /**    /**
477       \brief       \brief
478       Return a copy of the domain.       Return a copy of the domain.
# Line 452  class Data { Line 487  class Data {
487    */    */
488    ESCRIPT_DLL_API    ESCRIPT_DLL_API
489    inline    inline
490    int    unsigned int
491    getDataPointRank() const    getDataPointRank() const
492    {    {
493      return m_data->getPointDataView().getRank();      return m_data->getRank();
494    }    }
495    
496    /**    /**
# Line 493  class Data { Line 528  class Data {
528      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
529    }    }
530    
531    
532      /**
533        \brief
534        Return the number of values in the shape for this object.
535      */
536      ESCRIPT_DLL_API
537      int
538      getNoValues() const
539      {
540        return m_data->getNoValues();
541      }
542    
543    
544      /**
545         \brief
546         dumps the object into a netCDF file
547      */
548      ESCRIPT_DLL_API
549      void
550      dump(const std::string fileName) const;
551    
552     /**
553      \brief returns the values of the object as a list of tuples (one for each datapoint).
554    
555      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
556    If false, the result is a list of scalars [1, 2, ...]
557     */
558      ESCRIPT_DLL_API
559      const boost::python::object
560      toListOfTuples(bool scalarastuple=true);
561    
562    
563     /**
564        \brief
565        Return the sample data for the given sample no. This is not the
566        preferred interface but is provided for use by C code.
567        The bufferg parameter is only required for LazyData.
568        \param sampleNo - Input - the given sample no.
569        \return pointer to the sample data.
570    */
571      ESCRIPT_DLL_API
572      inline
573      const DataAbstract::ValueType::value_type*
574      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo);
575    
576    
577    /**    /**
578       \brief       \brief
579       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
580       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
581       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
582         \return pointer to the sample data.
583    */    */
584    ESCRIPT_DLL_API    ESCRIPT_DLL_API
585    inline    inline
586    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
587    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
588    {  
     return m_data->getSampleData(sampleNo);  
   }  
589    
590    /**    /**
591       \brief       \brief
# Line 523  class Data { Line 603  class Data {
603    
604    /**    /**
605       \brief       \brief
606       Assign the given value to the data-points referenced by the given       Return a reference into the DataVector which points to the specified data point.
607       reference number.       \param sampleNo - Input -
608         \param dataPointNo - Input -
      The value supplied is a python numarray object.  The data from this numarray  
      is unpacked into a DataArray, and this is used to set the corresponding  
      data-points in the underlying Data object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Input - value to assign to data-points associated with  
                             the given reference number.  
609    */    */
610    ESCRIPT_DLL_API    ESCRIPT_DLL_API
611    void    DataTypes::ValueType::const_reference
612    setRefValue(int ref,    getDataPointRO(int sampleNo, int dataPointNo);
               const boost::python::numeric::array& value);  
613    
614    /**    /**
615       \brief       \brief
616       Return the values associated with the data-points referenced by the given       Return a reference into the DataVector which points to the specified data point.
617       reference number.       \param sampleNo - Input -
618         \param dataPointNo - Input -
      The value supplied is a python numarray object. The data from the corresponding  
      data-points in this Data object are packed into the given numarray object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Output - object to receive values from data-points  
                              associated with the given reference number.  
619    */    */
620    ESCRIPT_DLL_API    ESCRIPT_DLL_API
621    void    DataTypes::ValueType::reference
622    getRefValue(int ref,    getDataPointRW(int sampleNo, int dataPointNo);
623                boost::python::numeric::array& value);  
624    
625    
626    /**    /**
627       \brief       \brief
628       Return a view into the data for the data point specified.       Return the offset for the given sample and point within the sample
      NOTE: Construction of the DataArrayView is a relatively expensive  
      operation.  
      \param sampleNo - Input -  
      \param dataPointNo - Input -  
629    */    */
630    ESCRIPT_DLL_API    ESCRIPT_DLL_API
631    inline    inline
632    DataArrayView    DataTypes::ValueType::size_type
633    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
634                 int dataPointNo)                 int dataPointNo)
635    {    {
636          return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
637    }    }
638    
639    /**    /**
# Line 584  class Data { Line 641  class Data {
641       Return a reference to the data point shape.       Return a reference to the data point shape.
642    */    */
643    ESCRIPT_DLL_API    ESCRIPT_DLL_API
644    const DataArrayView::ShapeType&    inline
645    getDataPointShape() const;    const DataTypes::ShapeType&
646      getDataPointShape() const
647      {
648        return m_data->getShape();
649      }
650    
651    /**    /**
652       \brief       \brief
# Line 609  class Data { Line 670  class Data {
670       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
671    */    */
672    ESCRIPT_DLL_API    ESCRIPT_DLL_API
673    DataArrayView::ValueType::size_type    DataTypes::ValueType::size_type
674    getLength() const;    getLength() const;
675    
676    /**    /**
677      \brief Return true if this object contains no samples.
678      This is not the same as isEmpty()
679      */
680      ESCRIPT_DLL_API
681      bool
682      hasNoSamples() const
683      {
684        return getLength()==0;
685      }
686    
687      /**
688       \brief       \brief
689       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
690       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
691       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
692         \param name - Input - name of tag.
693         \param value - Input - Value to associate with given key.
694      */
695      ESCRIPT_DLL_API
696      void
697      setTaggedValueByName(std::string name,
698                           const boost::python::object& value);
699    
700      /**
701         \brief
702         Assign the given value to the tag. Implicitly converts this
703         object to type DataTagged if it is constant.
704    
705       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
706       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
707      ==>*      ==>*
# Line 629  class Data { Line 714  class Data {
714    /**    /**
715       \brief       \brief
716       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
717       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
718       cannot be converted to a DataTagged object.  
719       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
720         \param pointshape - Input - The shape of the value parameter
721       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
722      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
723    */    */
724    ESCRIPT_DLL_API    ESCRIPT_DLL_API
725    void    void
726    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
727                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
728                            const DataTypes::ValueType& value,
729                int dataOffset=0);
730    
731    
732    
733    /**    /**
734      \brief      \brief
# Line 655  class Data { Line 745  class Data {
745    
746    /**    /**
747       \brief       \brief
748         set all values to zero
749         *
750      */
751      ESCRIPT_DLL_API
752      void
753      setToZero();
754    
755      /**
756         \brief
757       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
758       the result as a Data object.       the result as a Data object.
759       *       *
# Line 663  class Data { Line 762  class Data {
762    Data    Data
763    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
764    
765    
766      ESCRIPT_DLL_API
767      Data
768      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
769                           double undef, Data& B, double Bmin, double Bstep,bool check_boundaries);
770    
771      ESCRIPT_DLL_API
772      Data
773      interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
774                           double undef,bool check_boundaries);
775    
776    
777    
778    
779      ESCRIPT_DLL_API
780      Data
781      interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
782                            Data& B, double Bmin, double Bstep, double undef,bool check_boundaries);
783    
784      ESCRIPT_DLL_API
785      Data
786      interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
787                            double undef,bool check_boundaries);
788    
789    /**    /**
790       \brief       \brief
791       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 678  class Data { Line 801  class Data {
801    grad() const;    grad() const;
802    
803    /**    /**
804       \brief      \brief
805       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
806       *    */
807      ESCRIPT_DLL_API
808      boost::python::object
809      integrateToTuple_const() const;
810    
811    
812      /**
813        \brief
814         Calculate the integral over the function space domain as a python tuple.
815    */    */
816    ESCRIPT_DLL_API    ESCRIPT_DLL_API
817    boost::python::numeric::array    boost::python::object
818    integrate() const;    integrateToTuple();
819    
820    
821    
822    /**    /**
823       \brief       \brief
# Line 751  class Data { Line 884  class Data {
884    /**    /**
885       \brief       \brief
886       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
887       *  
888         The method is not const because lazy data needs to be expanded before Lsup can be computed.
889         The _const form can be used when the Data object is const, however this will only work for
890         Data which is not Lazy.
891    
892         For Data which contain no samples (or tagged Data for which no tags in use have a value)
893         zero is returned.
894    */    */
895    ESCRIPT_DLL_API    ESCRIPT_DLL_API
896    double    double
897    Lsup() const;    Lsup();
898    
   /**  
      \brief  
      Return the minimum absolute value of this Data object.  
      *  
   */  
899    ESCRIPT_DLL_API    ESCRIPT_DLL_API
900    double    double
901    Linf() const;    Lsup_const() const;
902    
903    
904    /**    /**
905       \brief       \brief
906       Return the maximum value of this Data object.       Return the maximum value of this Data object.
907       *  
908         The method is not const because lazy data needs to be expanded before sup can be computed.
909         The _const form can be used when the Data object is const, however this will only work for
910         Data which is not Lazy.
911    
912         For Data which contain no samples (or tagged Data for which no tags in use have a value)
913         a large negative value is returned.
914    */    */
915    ESCRIPT_DLL_API    ESCRIPT_DLL_API
916    double    double
917    sup() const;    sup();
918    
919      ESCRIPT_DLL_API
920      double
921      sup_const() const;
922    
923    
924    /**    /**
925       \brief       \brief
926       Return the minimum value of this Data object.       Return the minimum value of this Data object.
927       *  
928         The method is not const because lazy data needs to be expanded before inf can be computed.
929         The _const form can be used when the Data object is const, however this will only work for
930         Data which is not Lazy.
931    
932         For Data which contain no samples (or tagged Data for which no tags in use have a value)
933         a large positive value is returned.
934    */    */
935    ESCRIPT_DLL_API    ESCRIPT_DLL_API
936    double    double
937    inf() const;    inf();
938    
939      ESCRIPT_DLL_API
940      double
941      inf_const() const;
942    
943    
944    
945    /**    /**
946       \brief       \brief
# Line 814  class Data { Line 972  class Data {
972    /**    /**
973       \brief       \brief
974       Return the (sample number, data-point number) of the data point with       Return the (sample number, data-point number) of the data point with
975       the minimum value in this Data object.       the minimum component value in this Data object.
976         \note If you are working in python, please consider using Locator
977    instead of manually manipulating process and point IDs.
978    */    */
979    ESCRIPT_DLL_API    ESCRIPT_DLL_API
980    const boost::python::tuple    const boost::python::tuple
981    minGlobalDataPoint() const;    minGlobalDataPoint() const;
982    
983      /**
984         \brief
985         Return the (sample number, data-point number) of the data point with
986         the minimum component value in this Data object.
987         \note If you are working in python, please consider using Locator
988    instead of manually manipulating process and point IDs.
989      */
990    ESCRIPT_DLL_API    ESCRIPT_DLL_API
991    void    const boost::python::tuple
992    calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;    maxGlobalDataPoint() const;
993    
994    
995    
996    /**    /**
997       \brief       \brief
998       Return the sign of each data point of this Data object.       Return the sign of each data point of this Data object.
# Line 882  class Data { Line 1052  class Data {
1052    /**    /**
1053       \brief       \brief
1054       Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.       Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1055       the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than       the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1056       tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the       tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1057       first non-zero entry is positive.       first non-zero entry is positive.
1058       Currently this function is restricted to rank 2, square shape, and dimension 3       Currently this function is restricted to rank 2, square shape, and dimension 3
1059       *       *
# Line 1087  class Data { Line 1257  class Data {
1257    /**    /**
1258       \brief       \brief
1259       Return the given power of each data point of this boost python object.       Return the given power of each data point of this boost python object.
1260        
1261       \param right Input - the power to raise the object to.       \param right Input - the power to raise the object to.
1262       *       *
1263     */     */
# Line 1098  class Data { Line 1268  class Data {
1268    /**    /**
1269       \brief       \brief
1270       Return the given power of each data point of this boost python object.       Return the given power of each data point of this boost python object.
1271        
1272       \param left Input - the bases       \param left Input - the bases
1273       *       *
1274     */     */
# Line 1123  class Data { Line 1293  class Data {
1293    void    void
1294    saveVTK(std::string fileName) const;    saveVTK(std::string fileName) const;
1295    
1296    
1297    
1298    /**    /**
1299       \brief       \brief
1300       Overloaded operator +=       Overloaded operator +=
# Line 1134  class Data { Line 1306  class Data {
1306    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1307    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1308    
1309      ESCRIPT_DLL_API
1310      Data& operator=(const Data& other);
1311    
1312    /**    /**
1313       \brief       \brief
1314       Overloaded operator -=       Overloaded operator -=
# Line 1168  class Data { Line 1343  class Data {
1343    Data& operator/=(const boost::python::object& right);    Data& operator/=(const boost::python::object& right);
1344    
1345    /**    /**
1346        \brief return inverse of matricies.
1347      */
1348      ESCRIPT_DLL_API
1349      Data
1350      matrixInverse() const;
1351    
1352      /**
1353       \brief       \brief
1354       Returns true if this can be interpolated to functionspace.       Returns true if this can be interpolated to functionspace.
1355    */    */
# Line 1226  class Data { Line 1408  class Data {
1408    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1409    inline    inline
1410    void    void
1411    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1412    
1413    /**    /**
1414       \brief       \brief
# Line 1237  class Data { Line 1419  class Data {
1419    */    */
1420    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1421    Data    Data
1422    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1423    
1424    /**    /**
1425       \brief       \brief
# Line 1250  class Data { Line 1432  class Data {
1432    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1433    void    void
1434    setSlice(const Data& value,    setSlice(const Data& value,
1435             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1436    
1437    /**    /**
1438       \brief       \brief
1439       Archive the current Data object to the given file.       print the data values to stdout. Used for debugging
      \param fileName - Input - file to archive to.  
1440    */    */
1441    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1442    void    void
1443    archiveData(const std::string fileName);          print(void);
1444    
1445    /**    /**
1446       \brief       \brief
1447       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1448       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1449       Note - the current object must be of type DataEmpty.                   is returned
      \param fileName - Input - file to extract from.  
      \param fspace - Input - a suitable FunctionSpace descibing the data.  
1450    */    */
1451    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1452    void          int
1453    extractData(const std::string fileName,          get_MPIRank(void) const;
               const FunctionSpace& fspace);  
   
1454    
1455    /**    /**
1456       \brief       \brief
1457       print the data values to stdout. Used for debugging       return the MPI rank number of the local data
1458                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1459                     is returned
1460    */    */
1461    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1462    void          int
1463      print(void);          get_MPISize(void) const;
1464    
1465    /**    /**
1466       \brief       \brief
1467       return the MPI rank number of the local data       return the MPI rank number of the local data
1468           MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()                   MPI_COMM_WORLD is assumed and returned.
          is returned  
1469    */    */
1470    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1471      int          MPI_Comm
1472      get_MPIRank(void) const;          get_MPIComm(void) const;
1473    
1474    /**    /**
1475       \brief       \brief
1476       return the MPI rank number of the local data       return the object produced by the factory, which is a DataConstant or DataExpanded
1477           MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()      TODO Ownership of this object should be explained in doco.
          is returned  
1478    */    */
1479    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1480      int          DataAbstract*
1481      get_MPISize(void) const;          borrowData(void) const;
1482    
1483      ESCRIPT_DLL_API
1484            DataAbstract_ptr
1485            borrowDataPtr(void) const;
1486    
   /**  
      \brief  
      return the MPI rank number of the local data  
          MPI_COMM_WORLD is assumed and returned.  
   */  
1487    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1488      MPI_Comm          DataReady_ptr
1489      get_MPIComm(void) const;          borrowReadyPtr(void) const;
1490    
1491    
1492    
1493    /**    /**
1494       \brief       \brief
1495       return the object produced by the factory, which is a DataConstant or DataExpanded       Return a pointer to the beginning of the datapoint at the specified offset.
1496         TODO Eventually these should be inlined.
1497         \param i - position(offset) in the underlying datastructure
1498    */    */
1499    
1500      ESCRIPT_DLL_API
1501            DataTypes::ValueType::const_reference
1502            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1503    
1504    
1505    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1506      DataAbstract*          DataTypes::ValueType::reference
1507      borrowData(void) const;          getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1508    
1509    
1510    
1511   protected:   protected:
1512    
1513   private:   private:
1514    
1515    template <class BinaryOp>
1516      double
1517    #ifdef PASO_MPI
1518      lazyAlgWorker(double init, MPI_Op mpiop_type);
1519    #else
1520      lazyAlgWorker(double init);
1521    #endif
1522    
1523      double
1524      LsupWorker() const;
1525    
1526      double
1527      supWorker() const;
1528    
1529      double
1530      infWorker() const;
1531    
1532      boost::python::object
1533      integrateWorker() const;
1534    
1535      void
1536      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1537    
1538      void
1539      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1540    
1541      // For internal use in Data.cpp only!
1542      // other uses should call the main entry points and allow laziness
1543      Data
1544      minval_nonlazy() const;
1545    
1546      // For internal use in Data.cpp only!
1547      Data
1548      maxval_nonlazy() const;
1549    
1550    
1551    /**    /**
1552       \brief       \brief
1553       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1395  class Data { Line 1619  class Data {
1619       \brief       \brief
1620       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1621    */    */
1622    template <class IValueType>  
1623    void    void
1624    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1625             const DataTypes::ShapeType& shape,
1626               const FunctionSpace& what,               const FunctionSpace& what,
1627               bool expanded);               bool expanded);
1628    
1629      void
1630      initialise(const WrappedArray& value,
1631                     const FunctionSpace& what,
1632                     bool expanded);
1633    
1634    //    //
1635    // flag to protect the data object against any update    // flag to protect the data object against any update
1636    bool m_protected;    bool m_protected;
1637      mutable bool m_shared;
1638      bool m_lazy;
1639    
1640    //    //
1641    // pointer to the actual data object    // pointer to the actual data object
1642    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1643      DataAbstract_ptr m_data;
1644    
1645    // If possible please use getReadyPtr instead.
1646    // But see warning below.
1647      const DataReady*
1648      getReady() const;
1649    
1650      DataReady*
1651      getReady();
1652    
1653    
1654    // Be wary of using this for local operations since it (temporarily) increases reference count.
1655    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1656    // getReady() instead
1657      DataReady_ptr
1658      getReadyPtr();
1659    
1660      const_DataReady_ptr
1661      getReadyPtr() const;
1662    
1663    
1664      /**
1665       \brief Update the Data's shared flag
1666       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1667       For internal use only.
1668      */
1669      void updateShareStatus(bool nowshared) const
1670      {
1671        m_shared=nowshared;     // m_shared is mutable
1672      }
1673    
1674      // In the isShared() method below:
1675      // A problem would occur if m_data (the address pointed to) were being modified
1676      // while the call m_data->is_shared is being executed.
1677      //
1678      // Q: So why do I think this code can be thread safe/correct?
1679      // A: We need to make some assumptions.
1680      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1681      // 2. We assume that no constructions or assignments which will share previously unshared
1682      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1683    //    //
1684    // pointer to the internal profiling data    // This means that the only transition we need to consider, is when a previously shared object is
1685    struct profDataEntry *profData;    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1686      // In those cases the m_shared flag changes to false after m_data has completed changing.
1687      // For any threads executing before the flag switches they will assume the object is still shared.
1688      bool isShared() const
1689      {
1690        return m_shared;
1691    /*  if (m_shared) return true;
1692        if (m_data->isShared())        
1693        {                  
1694            updateShareStatus(true);
1695            return true;
1696        }
1697        return false;*/
1698      }
1699    
1700      void forceResolve()
1701      {
1702        if (isLazy())
1703        {
1704            #ifdef _OPENMP
1705            if (omp_in_parallel())
1706            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1707            throw DataException("Please do not call forceResolve() in a parallel region.");
1708            }
1709            #endif
1710            resolve();
1711        }
1712      }
1713    
1714      /**
1715      \brief if another object is sharing out member data make a copy to work with instead.
1716      This code should only be called from single threaded sections of code.
1717      */
1718      void exclusiveWrite()
1719      {
1720    #ifdef _OPENMP
1721      if (omp_in_parallel())
1722      {
1723    // *((int*)0)=17;
1724        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1725      }
1726    #endif
1727        forceResolve();
1728        if (isShared())
1729        {
1730            DataAbstract* t=m_data->deepCopy();
1731            set_m_data(DataAbstract_ptr(t));
1732        }
1733      }
1734    
1735      /**
1736      \brief checks if caller can have exclusive write to the object
1737      */
1738      void checkExclusiveWrite()
1739      {
1740        if  (isLazy() || isShared())
1741        {
1742            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1743        }
1744      }
1745    
1746      /**
1747      \brief Modify the data abstract hosted by this Data object
1748      For internal use only.
1749      Passing a pointer to null is permitted (do this in the destructor)
1750      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1751      */
1752      void set_m_data(DataAbstract_ptr p);
1753    
1754      friend class DataAbstract;        // To allow calls to updateShareStatus
1755    
1756  };  };
1757    
1758  template <class IValueType>  }   // end namespace escript
1759  void  
1760  Data::initialise(const IValueType& value,  
1761                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1762                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1763    // so that I can dynamic cast between them below.
1764    #include "DataReady.h"
1765    #include "DataLazy.h"
1766    
1767    namespace escript
1768  {  {
1769    //  
1770    // Construct a Data object of the appropriate type.  inline
1771    // Construct the object first as there seems to be a bug which causes  const DataReady*
1772    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1773    // within the shared_ptr constructor.  {
1774    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1775      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1776      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
1777      m_data=temp_data;  }
1778    } else {  
1779      DataAbstract* temp=new DataConstant(value,what);  inline
1780      boost::shared_ptr<DataAbstract> temp_data(temp);  DataReady*
1781      m_data=temp_data;  Data::getReady()
1782    }  {
1783       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1784       EsysAssert((dr!=0), "Error - casting to DataReady.");
1785       return dr;
1786    }
1787    
1788    // Be wary of using this for local operations since it (temporarily) increases reference count.
1789    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1790    // getReady() instead
1791    inline
1792    DataReady_ptr
1793    Data::getReadyPtr()
1794    {
1795       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1796       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1797       return dr;
1798  }  }
1799    
1800    
1801    inline
1802    const_DataReady_ptr
1803    Data::getReadyPtr() const
1804    {
1805       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1806       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1807       return dr;
1808    }
1809    
1810    inline
1811    DataAbstract::ValueType::value_type*
1812    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1813    {
1814       if (isLazy())
1815       {
1816        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1817       }
1818       return getReady()->getSampleDataRW(sampleNo);
1819    }
1820    
1821    inline
1822    const DataAbstract::ValueType::value_type*
1823    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo)
1824    {
1825       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1826       if (l!=0)
1827       {
1828        size_t offset=0;
1829        const DataTypes::ValueType* res=l->resolveSample(sampleNo,offset);
1830        return &((*res)[offset]);
1831       }
1832       return getReady()->getSampleDataRO(sampleNo);
1833    }
1834    
1835    
1836    
1837    /**
1838       Modify a filename for MPI parallel output to multiple files
1839    */
1840    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1841    
1842  /**  /**
1843     Binary Data object operators.     Binary Data object operators.
1844  */  */
1845  inline double rpow(double x,double y)  inline double rpow(double x,double y)
1846  {  {
1847      return pow(y,x);      return pow(y,x);
1848  };  }
1849    
1850  /**  /**
1851    \brief    \brief
# Line 1537  ESCRIPT_DLL_API Data operator*(const boo Line 1939  ESCRIPT_DLL_API Data operator*(const boo
1939  */  */
1940  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1941    
1942    
1943    
1944  /**  /**
1945    \brief    \brief
1946    Output operator    Output operator
# Line 1546  ESCRIPT_DLL_API std::ostream& operator<< Line 1950  ESCRIPT_DLL_API std::ostream& operator<<
1950  /**  /**
1951    \brief    \brief
1952    Compute a tensor product of two Data objects    Compute a tensor product of two Data objects
1953    \param arg0 - Input - Data object    \param arg_0 - Input - Data object
1954    \param arg1 - Input - Data object    \param arg_1 - Input - Data object
1955    \param axis_offset - Input - axis offset    \param axis_offset - Input - axis offset
1956    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1    \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1957  */  */
1958  ESCRIPT_DLL_API  ESCRIPT_DLL_API
1959  Data  Data
1960  C_GeneralTensorProduct(Data& arg0,  C_GeneralTensorProduct(Data& arg_0,
1961                       Data& arg1,                       Data& arg_1,
1962                       int axis_offset=0,                       int axis_offset=0,
1963                       int transpose=0);                       int transpose=0);
1964    
1965  /**  /**
1966    \brief    \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);  
   
 /**  
   \brief  
1967    Perform the given binary operation with this and right as operands.    Perform the given binary operation with this and right as operands.
1968    Right is a Data object.    Right is a Data object.
1969  */  */
# Line 1579  Data::binaryOp(const Data& right, Line 1975  Data::binaryOp(const Data& right,
1975  {  {
1976     //     //
1977     // 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
1978     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1979       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.");
1980     }     }
1981    
1982       if (isLazy() || right.isLazy())
1983       {
1984         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1985       }
1986     //     //
1987     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1988     Data tempRight(right);     Data tempRight(right);
1989    
1990     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1991       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1992         //         //
1993         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1994         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1995       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1996         //         //
1997         // interpolate onto the RHS function space         // interpolate onto the RHS function space
1998         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
1999         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2000           set_m_data(tempLeft.m_data);
2001       }       }
2002     }     }
2003     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1610  Data::binaryOp(const Data& right, Line 2013  Data::binaryOp(const Data& right,
2013       // of any data type       // of any data type
2014       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2015       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2016       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2017     } else if (isTagged()) {     } else if (isTagged()) {
2018       //       //
2019       // 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 1632  Data::binaryOp(const Data& right, Line 2035  Data::binaryOp(const Data& right,
2035       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");       EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2036       escript::binaryOp(*leftC,*rightC,operation);       escript::binaryOp(*leftC,*rightC,operation);
2037     }     }
    #if defined DOPROF  
    profData->binary++;  
    #endif  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on other and return the result.  
   Given operation is performed on each element of each data point, thus  
   argument object is a rank n Data object, and returned object is a rank n  
   Data object.  
   Calls Data::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 Data  
 unaryOp(const Data& other,  
         UnaryFunction operation)  
 {  
   Data result;  
   result.copy(other);  
   result.unaryOp(operation);  
   return result;  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on this.  
   Given operation is performed on each element of each data point.  
   Calls escript::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 void  
 Data::unaryOp(UnaryFunction operation)  
 {  
   if (isExpanded()) {  
     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isTagged()) {  
     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isConstant()) {  
     DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");  
     escript::unaryOp(*leftC,operation);  
   }  
2038  }  }
2039    
2040  /**  /**
# Line 1707  Data::algorithm(BinaryFunction operation Line 2061  Data::algorithm(BinaryFunction operation
2061      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2062      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2063      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
2064      } else if (isEmpty()) {
2065        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2066      } else if (isLazy()) {
2067        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2068      } else {
2069        throw DataException("Error - Data encapsulates an unknown type.");
2070    }    }
   return 0;  
2071  }  }
2072    
2073  /**  /**
2074    \brief    \brief
2075    Perform the given data point reduction algorithm on data and return the result.    Perform the given data point reduction algorithm on data and return the result.
2076    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2077    thus argument object is a rank n Data object, and returned object is a    thus argument object is a rank n Data object, and returned object is a
2078    rank 0 Data object.    rank 0 Data object.
2079    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2080  */  */
# Line 1724  inline Line 2083  inline
2083  Data  Data
2084  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2085  {  {
2086    if (isExpanded()) {    if (isEmpty()) {
2087      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2088      }
2089      else if (isExpanded()) {
2090        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2091      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2092      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2093      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2094      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2095      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2096      return result;      return result;
2097    } else if (isTagged()) {    }
2098      else if (isTagged()) {
2099      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());  
2100      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2101      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2102        defval[0]=0;
2103        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2104      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2105      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
2106    } else if (isConstant()) {    }
2107      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
2108        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2109      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2110      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2111      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2112      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2113      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2114      return result;      return result;
2115      } else if (isLazy()) {
2116        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2117      } else {
2118        throw DataException("Error - Data encapsulates an unknown type.");
2119      }
2120    }
2121    
2122    /**
2123      \brief
2124      Compute a tensor operation with two Data objects
2125      \param arg_0 - Input - Data object
2126      \param arg_1 - Input - Data object
2127      \param operation - Input - Binary op functor
2128    */
2129    template <typename BinaryFunction>
2130    inline
2131    Data
2132    C_TensorBinaryOperation(Data const &arg_0,
2133                            Data const &arg_1,
2134                            BinaryFunction operation)
2135    {
2136      if (arg_0.isEmpty() || arg_1.isEmpty())
2137      {
2138         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2139      }
2140      if (arg_0.isLazy() || arg_1.isLazy())
2141      {
2142         throw DataException("Error - Operations not permitted on lazy data.");
2143      }
2144      // Interpolate if necessary and find an appropriate function space
2145      Data arg_0_Z, arg_1_Z;
2146      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2147        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2148          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2149          arg_1_Z = Data(arg_1);
2150        }
2151        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2152          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2153          arg_0_Z =Data(arg_0);
2154        }
2155        else {
2156          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2157        }
2158      } else {
2159          arg_0_Z = Data(arg_0);
2160          arg_1_Z = Data(arg_1);
2161      }
2162      // Get rank and shape of inputs
2163      int rank0 = arg_0_Z.getDataPointRank();
2164      int rank1 = arg_1_Z.getDataPointRank();
2165      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2166      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2167      int size0 = arg_0_Z.getDataPointSize();
2168      int size1 = arg_1_Z.getDataPointSize();
2169      // Declare output Data object
2170      Data res;
2171    
2172      if (shape0 == shape1) {
2173        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2174          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2175          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2176          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2177          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2178    
2179          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2180        }
2181        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2182    
2183          // Prepare the DataConstant input
2184          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2185    
2186          // Borrow DataTagged input from Data object
2187          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2188    
2189          // Prepare a DataTagged output 2
2190          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2191          res.tag();
2192          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2193    
2194          // Prepare offset into DataConstant
2195          int offset_0 = tmp_0->getPointOffset(0,0);
2196          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2197    
2198          // Get the pointers to the actual data
2199          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2200          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2201    
2202          // Compute a result for the default
2203          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2204          // Compute a result for each tag
2205          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2206          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2207          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2208            tmp_2->addTag(i->first);
2209            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2210            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2211    
2212            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2213          }
2214    
2215        }
2216        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2217          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2218          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2219          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2220          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2221    
2222          int sampleNo_1,dataPointNo_1;
2223          int numSamples_1 = arg_1_Z.getNumSamples();
2224          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2225          int offset_0 = tmp_0->getPointOffset(0,0);
2226          res.requireWrite();
2227          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2228          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2229            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2230              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2231              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2232              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2233              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2234              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2235              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2236            }
2237          }
2238    
2239        }
2240        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2241          // Borrow DataTagged input from Data object
2242          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2243    
2244          // Prepare the DataConstant input
2245          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2246    
2247          // Prepare a DataTagged output 2
2248          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2249          res.tag();
2250          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2251    
2252          // Prepare offset into DataConstant
2253          int offset_1 = tmp_1->getPointOffset(0,0);
2254    
2255          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2256          // Get the pointers to the actual data
2257          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2258          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2259          // Compute a result for the default
2260          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2261          // Compute a result for each tag
2262          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2263          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2264          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2265            tmp_2->addTag(i->first);
2266            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2267            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2268            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2269          }
2270    
2271        }
2272        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2273          // Borrow DataTagged input from Data object
2274          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2275    
2276          // Borrow DataTagged input from Data object
2277          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2278    
2279          // Prepare a DataTagged output 2
2280          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2281          res.tag();        // DataTagged output
2282          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2283    
2284          // Get the pointers to the actual data
2285          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2286          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2287          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2288    
2289          // Compute a result for the default
2290          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2291          // Merge the tags
2292          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2293          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2294          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2295          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2296            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2297          }
2298          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2299            tmp_2->addTag(i->first);
2300          }
2301          // Compute a result for each tag
2302          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2303          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2304    
2305            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2306            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2307            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2308    
2309            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2310          }
2311    
2312        }
2313        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2314          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2315          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2316          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2317          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2318          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2319    
2320          int sampleNo_0,dataPointNo_0;
2321          int numSamples_0 = arg_0_Z.getNumSamples();
2322          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2323          res.requireWrite();
2324          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2325          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2326            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2327            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2328            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2329              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2330              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2331              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2332              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2333              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2334            }
2335          }
2336    
2337        }
2338        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2339          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2340          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2341          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2342          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2343    
2344          int sampleNo_0,dataPointNo_0;
2345          int numSamples_0 = arg_0_Z.getNumSamples();
2346          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2347          int offset_1 = tmp_1->getPointOffset(0,0);
2348          res.requireWrite();
2349          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2350          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2351            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2352              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2353              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2354    
2355              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2356              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2357              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2358    
2359    
2360              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2361            }
2362          }
2363    
2364        }
2365        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2366          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2367          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2368          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2369          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2370          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2371    
2372          int sampleNo_0,dataPointNo_0;
2373          int numSamples_0 = arg_0_Z.getNumSamples();
2374          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2375          res.requireWrite();
2376          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2377          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2378            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2379            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2380            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2381              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2382              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2383              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2384              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2385              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2386            }
2387          }
2388    
2389        }
2390        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2391          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2392          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2393          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2394          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2395          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2396    
2397          int sampleNo_0,dataPointNo_0;
2398          int numSamples_0 = arg_0_Z.getNumSamples();
2399          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2400          res.requireWrite();
2401          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2402          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2403          dataPointNo_0=0;
2404    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2405              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2406              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2407              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2408              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2409              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2410              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2411              tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2412    //       }
2413          }
2414    
2415        }
2416        else {
2417          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2418        }
2419    
2420      } else if (0 == rank0) {
2421        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2422          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2423          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2424          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2425          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2426          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2427        }
2428        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2429    
2430          // Prepare the DataConstant input
2431          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2432    
2433          // Borrow DataTagged input from Data object
2434          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2435    
2436          // Prepare a DataTagged output 2
2437          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2438          res.tag();
2439          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2440    
2441          // Prepare offset into DataConstant
2442          int offset_0 = tmp_0->getPointOffset(0,0);
2443          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2444    
2445          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2446          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2447    
2448          // Compute a result for the default
2449          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2450          // Compute a result for each tag
2451          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2452          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2453          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2454            tmp_2->addTag(i->first);
2455            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2456            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2457            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2458          }
2459    
2460        }
2461        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2462    
2463          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2464          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2465          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2466          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2467    
2468          int sampleNo_1,dataPointNo_1;
2469          int numSamples_1 = arg_1_Z.getNumSamples();
2470          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2471          int offset_0 = tmp_0->getPointOffset(0,0);
2472          res.requireWrite();
2473          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2474          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2475            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2476              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2477              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2478              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2479              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2480              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2481              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2482    
2483            }
2484          }
2485    
2486        }
2487        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2488    
2489          // Borrow DataTagged input from Data object
2490          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2491    
2492          // Prepare the DataConstant input
2493          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2494    
2495          // Prepare a DataTagged output 2
2496          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2497          res.tag();
2498          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2499    
2500          // Prepare offset into DataConstant
2501          int offset_1 = tmp_1->getPointOffset(0,0);
2502          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2503    
2504          // Get the pointers to the actual data
2505          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2506          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2507    
2508    
2509          // Compute a result for the default
2510          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2511          // Compute a result for each tag
2512          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2513          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2514          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2515            tmp_2->addTag(i->first);
2516            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2517            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2518    
2519            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2520          }
2521    
2522        }
2523        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2524    
2525          // Borrow DataTagged input from Data object
2526          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2527    
2528          // Borrow DataTagged input from Data object
2529          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2530    
2531          // Prepare a DataTagged output 2
2532          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2533          res.tag();        // DataTagged output
2534          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2535    
2536          // Get the pointers to the actual data
2537          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2538          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2539          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2540    
2541          // Compute a result for the default
2542          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2543          // Merge the tags
2544          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2545          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2546          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2547          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2548            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2549          }
2550          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2551            tmp_2->addTag(i->first);
2552          }
2553          // Compute a result for each tag
2554          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2555          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2556            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2557            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2558            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2559    
2560            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2561          }
2562    
2563        }
2564        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2565    
2566          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2567          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2568          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2569          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2570          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2571    
2572          int sampleNo_0,dataPointNo_0;
2573          int numSamples_0 = arg_0_Z.getNumSamples();
2574          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2575          res.requireWrite();
2576          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2577          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2578            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2579            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2580            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2581              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2582              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2583              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2584              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2585              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2586            }
2587          }
2588    
2589        }
2590        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2591          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2592          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2593          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2594          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2595    
2596          int sampleNo_0,dataPointNo_0;
2597          int numSamples_0 = arg_0_Z.getNumSamples();
2598          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2599          int offset_1 = tmp_1->getPointOffset(0,0);
2600          res.requireWrite();
2601          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2602          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2603            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2604              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2605              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2606              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2607              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2608              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2609              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2610            }
2611          }
2612    
2613    
2614        }
2615        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2616    
2617          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2618          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2619          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2620          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2621          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2622    
2623          int sampleNo_0,dataPointNo_0;
2624          int numSamples_0 = arg_0_Z.getNumSamples();
2625          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2626          res.requireWrite();
2627          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2628          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2629            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2630            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2631            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2632              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2633              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2634              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2635              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2636              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2637            }
2638          }
2639    
2640        }
2641        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2642    
2643          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2644          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2645          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2646          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2647          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2648    
2649          int sampleNo_0,dataPointNo_0;
2650          int numSamples_0 = arg_0_Z.getNumSamples();
2651          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2652          res.requireWrite();
2653          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2654          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2655            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2656              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2657              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2658              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2659              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2660              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2661              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2662              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2663            }
2664          }
2665    
2666        }
2667        else {
2668          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2669        }
2670    
2671      } else if (0 == rank1) {
2672        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2673          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2674          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2675          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2676          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2677          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2678        }
2679        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2680    
2681          // Prepare the DataConstant input
2682          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2683    
2684          // Borrow DataTagged input from Data object
2685          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2686    
2687          // Prepare a DataTagged output 2
2688          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2689          res.tag();
2690          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2691    
2692          // Prepare offset into DataConstant
2693          int offset_0 = tmp_0->getPointOffset(0,0);
2694          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2695    
2696          //Get the pointers to the actual data
2697          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2698          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2699    
2700          // Compute a result for the default
2701          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2702          // Compute a result for each tag
2703          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2704          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2705          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2706            tmp_2->addTag(i->first);
2707            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2708            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2709            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2710          }
2711        }
2712        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2713    
2714          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2715          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2716          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2717          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2718    
2719          int sampleNo_1,dataPointNo_1;
2720          int numSamples_1 = arg_1_Z.getNumSamples();
2721          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2722          int offset_0 = tmp_0->getPointOffset(0,0);
2723          res.requireWrite();
2724          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2725          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2726            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2727              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2728              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2729              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2730              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2731              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2732              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2733            }
2734          }
2735    
2736        }
2737        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2738    
2739          // Borrow DataTagged input from Data object
2740          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2741    
2742          // Prepare the DataConstant input
2743          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2744    
2745          // Prepare a DataTagged output 2
2746          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2747          res.tag();
2748          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2749    
2750          // Prepare offset into DataConstant
2751          int offset_1 = tmp_1->getPointOffset(0,0);
2752          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2753          // Get the pointers to the actual data
2754          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2755          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2756          // Compute a result for the default
2757          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2758          // Compute a result for each tag
2759          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2760          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2761          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2762            tmp_2->addTag(i->first);
2763            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2764            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2765            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2766          }
2767    
2768        }
2769        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2770    
2771          // Borrow DataTagged input from Data object
2772          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2773    
2774          // Borrow DataTagged input from Data object
2775          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2776    
2777          // Prepare a DataTagged output 2
2778          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2779          res.tag();        // DataTagged output
2780          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2781    
2782          // Get the pointers to the actual data
2783          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2784          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2785          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2786    
2787          // Compute a result for the default
2788          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2789          // Merge the tags
2790          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2791          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2792          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2793          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2794            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2795          }
2796          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2797            tmp_2->addTag(i->first);
2798          }
2799          // Compute a result for each tag
2800          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2801          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2802            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2803            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2804            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2805            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2806          }
2807    
2808        }
2809        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2810    
2811          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2812          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2813          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2814          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2815          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2816    
2817          int sampleNo_0,dataPointNo_0;
2818          int numSamples_0 = arg_0_Z.getNumSamples();
2819          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2820          res.requireWrite();
2821          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2822          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2823            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2824            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2825            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2826              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2827              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2828              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2829              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2830              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2831            }
2832          }
2833    
2834        }
2835        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2836          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2837          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2838          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2839          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2840    
2841          int sampleNo_0,dataPointNo_0;
2842          int numSamples_0 = arg_0_Z.getNumSamples();
2843          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2844          int offset_1 = tmp_1->getPointOffset(0,0);
2845          res.requireWrite();
2846          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2847          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2848            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2849              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2850              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2851              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2852              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2853              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2854              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2855            }
2856          }
2857    
2858    
2859        }
2860        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2861    
2862          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2863          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2864          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2865          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2866          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2867    
2868          int sampleNo_0,dataPointNo_0;
2869          int numSamples_0 = arg_0_Z.getNumSamples();
2870          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2871          res.requireWrite();
2872          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2873          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2874            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2875            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2876            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2877              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2878              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2879              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2880              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2881              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2882            }
2883          }
2884    
2885        }
2886        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2887    
2888          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2889          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2890          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2891          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2892          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2893    
2894          int sampleNo_0,dataPointNo_0;
2895          int numSamples_0 = arg_0_Z.getNumSamples();
2896          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2897          res.requireWrite();
2898          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2899          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2900            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2901              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2902              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2903              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2904              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2905              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2906              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2907              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2908            }
2909          }
2910    
2911        }
2912        else {
2913          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2914        }
2915    
2916      } else {
2917        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2918      }
2919    
2920      return res;
2921    }
2922    
2923    template <typename UnaryFunction>
2924    Data
2925    C_TensorUnaryOperation(Data const &arg_0,
2926                           UnaryFunction operation)
2927    {
2928      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2929      {
2930         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2931      }
2932      if (arg_0.isLazy())
2933      {
2934         throw DataException("Error - Operations not permitted on lazy data.");
2935      }
2936      // Interpolate if necessary and find an appropriate function space
2937      Data arg_0_Z = Data(arg_0);
2938    
2939      // Get rank and shape of inputs
2940      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2941      int size0 = arg_0_Z.getDataPointSize();
2942    
2943      // Declare output Data object
2944      Data res;
2945    
2946      if (arg_0_Z.isConstant()) {
2947        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2948        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2949        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2950        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2951      }
2952      else if (arg_0_Z.isTagged()) {
2953    
2954        // Borrow DataTagged input from Data object
2955        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2956    
2957        // Prepare a DataTagged output 2
2958        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2959        res.tag();
2960        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2961    
2962        // Get the pointers to the actual data
2963        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2964        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2965        // Compute a result for the default
2966        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2967        // Compute a result for each tag
2968        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2969        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2970        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2971          tmp_2->addTag(i->first);
2972          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2973          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2974          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2975        }
2976    
2977    }    }
2978    Data falseRetVal; // to keep compiler quiet    else if (arg_0_Z.isExpanded()) {
2979    return falseRetVal;  
2980        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2981        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2982        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2983    
2984        int sampleNo_0,dataPointNo_0;
2985        int numSamples_0 = arg_0_Z.getNumSamples();
2986        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2987        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2988        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2989        dataPointNo_0=0;
2990    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2991            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2992            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2993            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2994            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2995            tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
2996    //      }
2997        }
2998      }
2999      else {
3000        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3001      }
3002    
3003      return res;
3004  }  }
3005    
3006  }  }

Legend:
Removed from v.922  
changed lines
  Added in v.2770

  ViewVC Help
Powered by ViewVC 1.1.26