Data.h

Go to the documentation of this file.

/*******************************************************
*
* Copyright (c) 2003-2010 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


#ifndef DATA_H
#define DATA_H
#include "system_dep.h"

#include "DataTypes.h"
#include "DataAbstract.h"
#include "DataAlgorithm.h"
#include "FunctionSpace.h"
#include "BinaryOp.h"
#include "UnaryOp.h"
#include "DataException.h"



extern "C" {
#include "DataC.h"
//#include <omp.h>
}

#ifdef _OPENMP
#include <omp.h>
#endif

#include "esysUtils/Esys_MPI.h"
#include <string>
#include <algorithm>
#include <sstream>

#include <boost/shared_ptr.hpp>
#include <boost/python/object.hpp>
#include <boost/python/tuple.hpp>

namespace escript {

//
// Forward declaration for various implementations of Data.
class DataConstant;
class DataTagged;
class DataExpanded;
class DataLazy;

class Data {

  public:

  // These typedefs allow function names to be cast to pointers
  // to functions of the appropriate type when calling unaryOp etc.
  typedef double (*UnaryDFunPtr)(double);
  typedef double (*BinaryDFunPtr)(double,double);


  ESCRIPT_DLL_API
  Data();

  ESCRIPT_DLL_API
  Data(const Data& inData);

  ESCRIPT_DLL_API
  Data(const Data& inData,
       const FunctionSpace& what);

  ESCRIPT_DLL_API
  Data(const DataTypes::ValueType& value,
         const DataTypes::ShapeType& shape,
                 const FunctionSpace& what=FunctionSpace(),
                 bool expanded=false);

  ESCRIPT_DLL_API
  Data(double value,
       const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
       const FunctionSpace& what=FunctionSpace(),
       bool expanded=false);

  ESCRIPT_DLL_API
  Data(const Data& inData,
       const DataTypes::RegionType& region);

  ESCRIPT_DLL_API
  Data(const boost::python::object& value,
       const FunctionSpace& what=FunctionSpace(),
       bool expanded=false);

  ESCRIPT_DLL_API     
  Data(const WrappedArray& w, const FunctionSpace& what,
           bool expanded=false);       
       

  ESCRIPT_DLL_API
  Data(const boost::python::object& value,
       const Data& other);

  ESCRIPT_DLL_API
  Data(double value,
       const boost::python::tuple& shape=boost::python::make_tuple(),
       const FunctionSpace& what=FunctionSpace(),
       bool expanded=false);

  ESCRIPT_DLL_API
  explicit Data(DataAbstract* underlyingdata);

  ESCRIPT_DLL_API
  explicit Data(DataAbstract_ptr underlyingdata);

  ESCRIPT_DLL_API
  ~Data();

  ESCRIPT_DLL_API
  void
  copy(const Data& other);

  ESCRIPT_DLL_API
  Data
  copySelf();


  ESCRIPT_DLL_API
  Data
  delay();

  ESCRIPT_DLL_API
  void 
  delaySelf();


  ESCRIPT_DLL_API
  void
  setProtection();

  ESCRIPT_DLL_API
  bool
  isProtected() const;


  ESCRIPT_DLL_API
  const boost::python::object
  getValueOfDataPointAsTuple(int dataPointNo);

  ESCRIPT_DLL_API
  void
  setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);

  ESCRIPT_DLL_API
  void
  setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);

  ESCRIPT_DLL_API
  void
  setValueOfDataPoint(int dataPointNo, const double);

  ESCRIPT_DLL_API
  const boost::python::object
  getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);

  
  ESCRIPT_DLL_API
  void
  setTupleForGlobalDataPoint(int id, int proc, boost::python::object);
  
  ESCRIPT_DLL_API
  int
  getTagNumber(int dpno);

  ESCRIPT_DLL_API
  escriptDataC
  getDataC();



  ESCRIPT_DLL_API
  escriptDataC
  getDataC() const;


  ESCRIPT_DLL_API
  std::string
  toString() const;

  ESCRIPT_DLL_API
  void
  expand();

  ESCRIPT_DLL_API
  void
  tag();

  ESCRIPT_DLL_API
  void
  resolve();


  ESCRIPT_DLL_API
  void
  requireWrite();

  ESCRIPT_DLL_API
  bool
  isExpanded() const;

  ESCRIPT_DLL_API
  bool
  actsExpanded() const;
  

  ESCRIPT_DLL_API
  bool
  isTagged() const;

  ESCRIPT_DLL_API
  bool
  isConstant() const;

  ESCRIPT_DLL_API
  bool
  isLazy() const;

  ESCRIPT_DLL_API
  bool
  isReady() const;

  ESCRIPT_DLL_API
  bool
  isEmpty() const;

  ESCRIPT_DLL_API
  inline
  const FunctionSpace&
  getFunctionSpace() const
  {
    return m_data->getFunctionSpace();
  }

  ESCRIPT_DLL_API
  const FunctionSpace
  getCopyOfFunctionSpace() const;

  ESCRIPT_DLL_API
  inline
//   const AbstractDomain&
  const_Domain_ptr
  getDomain() const
  {
     return getFunctionSpace().getDomain();
  }


  ESCRIPT_DLL_API
  inline
//   const AbstractDomain&
  Domain_ptr
  getDomainPython() const
  {
     return getFunctionSpace().getDomainPython();
  }

  ESCRIPT_DLL_API
  const AbstractDomain
  getCopyOfDomain() const;

  ESCRIPT_DLL_API
  inline
  unsigned int
  getDataPointRank() const
  {
    return m_data->getRank();
  }

  ESCRIPT_DLL_API
  inline
  int
  getNumDataPoints() const
  {
    return getNumSamples() * getNumDataPointsPerSample();
  }
  ESCRIPT_DLL_API
  inline
  int
  getNumSamples() const
  {
    return m_data->getNumSamples();
  }

  ESCRIPT_DLL_API
  inline
  int
  getNumDataPointsPerSample() const
  {
    return m_data->getNumDPPSample();
  }


  ESCRIPT_DLL_API
  int
  getNoValues() const
  {
    return m_data->getNoValues();
  }


  ESCRIPT_DLL_API
  void
  dump(const std::string fileName) const;

  ESCRIPT_DLL_API
  const boost::python::object
  toListOfTuples(bool scalarastuple=true);


  ESCRIPT_DLL_API
  inline
  const DataAbstract::ValueType::value_type*
  getSampleDataRO(DataAbstract::ValueType::size_type sampleNo);


  ESCRIPT_DLL_API
  inline
  DataAbstract::ValueType::value_type*
  getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);


  ESCRIPT_DLL_API
  inline
  DataAbstract::ValueType::value_type*
  getSampleDataByTag(int tag)
  {
    return m_data->getSampleDataByTag(tag);
  }

  ESCRIPT_DLL_API
  DataTypes::ValueType::const_reference
  getDataPointRO(int sampleNo, int dataPointNo);

  ESCRIPT_DLL_API
  DataTypes::ValueType::reference
  getDataPointRW(int sampleNo, int dataPointNo);



  ESCRIPT_DLL_API
  inline
  DataTypes::ValueType::size_type
  getDataOffset(int sampleNo,
               int dataPointNo)
  {
      return m_data->getPointOffset(sampleNo,dataPointNo);
  }

  ESCRIPT_DLL_API
  inline
  const DataTypes::ShapeType&
  getDataPointShape() const
  {
    return m_data->getShape();
  }

  ESCRIPT_DLL_API
  const boost::python::tuple
  getShapeTuple() const;

  ESCRIPT_DLL_API
  int
  getDataPointSize() const;

  ESCRIPT_DLL_API
  DataTypes::ValueType::size_type
  getLength() const;

  ESCRIPT_DLL_API
  bool
  hasNoSamples() const
  {
    return getLength()==0;
  }

  ESCRIPT_DLL_API
  void
  setTaggedValueByName(std::string name,
                       const boost::python::object& value);

  ESCRIPT_DLL_API
  void
  setTaggedValue(int tagKey,
                 const boost::python::object& value);

  ESCRIPT_DLL_API
  void
  setTaggedValueFromCPP(int tagKey,
            const DataTypes::ShapeType& pointshape,
                        const DataTypes::ValueType& value,
            int dataOffset=0);



  ESCRIPT_DLL_API
  void
  copyWithMask(const Data& other,
               const Data& mask);

  ESCRIPT_DLL_API
  void
  setToZero();

  ESCRIPT_DLL_API
  Data
  interpolate(const FunctionSpace& functionspace) const;

  ESCRIPT_DLL_API
  Data
  interpolateFromTable3D(const WrappedArray& table, double Amin, double Astep,
                       double undef, Data& B, double Bmin, double Bstep, Data& C, 
            double Cmin, double Cstep, bool check_boundaries);

  ESCRIPT_DLL_API
  Data
  interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
                       double undef, Data& B, double Bmin, double Bstep,bool check_boundaries);

  ESCRIPT_DLL_API
  Data
  interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
                       double undef,bool check_boundaries);


  ESCRIPT_DLL_API
  Data
  interpolateFromTable3DP(boost::python::object table, double Amin, double Astep,
                        Data& B, double Bmin, double Bstep, Data& C, double Cmin, double Cstep, double undef,bool check_boundaries);


  ESCRIPT_DLL_API
  Data
  interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
                        Data& B, double Bmin, double Bstep, double undef,bool check_boundaries);

  ESCRIPT_DLL_API
  Data
  interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
                        double undef,bool check_boundaries);

  ESCRIPT_DLL_API
  Data
  gradOn(const FunctionSpace& functionspace) const;

  ESCRIPT_DLL_API
  Data
  grad() const;

  ESCRIPT_DLL_API
  boost::python::object
  integrateToTuple_const() const;


  ESCRIPT_DLL_API
  boost::python::object
  integrateToTuple();



  ESCRIPT_DLL_API
  Data
  oneOver() const;
  ESCRIPT_DLL_API
  Data
  wherePositive() const;

  ESCRIPT_DLL_API
  Data
  whereNegative() const;

  ESCRIPT_DLL_API
  Data
  whereNonNegative() const;

  ESCRIPT_DLL_API
  Data
  whereNonPositive() const;

  ESCRIPT_DLL_API
  Data
  whereZero(double tol=0.0) const;

  ESCRIPT_DLL_API
  Data
  whereNonZero(double tol=0.0) const;

  ESCRIPT_DLL_API
  double
  Lsup();

  ESCRIPT_DLL_API
  double
  Lsup_const() const;


  ESCRIPT_DLL_API
  double
  sup();

  ESCRIPT_DLL_API
  double
  sup_const() const;


  ESCRIPT_DLL_API
  double
  inf();

  ESCRIPT_DLL_API
  double
  inf_const() const;



  ESCRIPT_DLL_API
  Data
  abs() const;

  ESCRIPT_DLL_API
  Data
  maxval() const;

  ESCRIPT_DLL_API
  Data
  minval() const;

  ESCRIPT_DLL_API
  const boost::python::tuple
  minGlobalDataPoint() const;

  ESCRIPT_DLL_API
  const boost::python::tuple
  maxGlobalDataPoint() const;



  ESCRIPT_DLL_API
  Data
  sign() const;

  ESCRIPT_DLL_API
  Data
  symmetric() const;

  ESCRIPT_DLL_API
  Data
  nonsymmetric() const;

  ESCRIPT_DLL_API
  Data
  trace(int axis_offset) const;

  ESCRIPT_DLL_API
  Data
  transpose(int axis_offset) const;

  ESCRIPT_DLL_API
  Data
  eigenvalues() const;

  ESCRIPT_DLL_API
  const boost::python::tuple
  eigenvalues_and_eigenvectors(const double tol=1.e-12) const;

  ESCRIPT_DLL_API
  Data
  swapaxes(const int axis0, const int axis1) const;

  ESCRIPT_DLL_API
  Data
  erf() const;

  ESCRIPT_DLL_API
  Data
  sin() const;

  ESCRIPT_DLL_API
  Data
  cos() const;

  ESCRIPT_DLL_API
  Data
  tan() const;

  ESCRIPT_DLL_API
  Data
  asin() const;

  ESCRIPT_DLL_API
  Data
  acos() const;

  ESCRIPT_DLL_API
  Data
  atan() const;

  ESCRIPT_DLL_API
  Data
  sinh() const;

  ESCRIPT_DLL_API
  Data
  cosh() const;

  ESCRIPT_DLL_API
  Data
  tanh() const;

  ESCRIPT_DLL_API
  Data
  asinh() const;

  ESCRIPT_DLL_API
  Data
  acosh() const;

  ESCRIPT_DLL_API
  Data
  atanh() const;

  ESCRIPT_DLL_API
  Data
  log10() const;

  ESCRIPT_DLL_API
  Data
  log() const;

  ESCRIPT_DLL_API
  Data
  exp() const;

  ESCRIPT_DLL_API
  Data
  sqrt() const;

  ESCRIPT_DLL_API
  Data
  neg() const;

  ESCRIPT_DLL_API
  Data
  pos() const;

  ESCRIPT_DLL_API
  Data
  powD(const Data& right) const;

  ESCRIPT_DLL_API
  Data
  powO(const boost::python::object& right) const;

  ESCRIPT_DLL_API
  Data
  rpowO(const boost::python::object& left) const;

  ESCRIPT_DLL_API
  void
  saveDX(std::string fileName) const;

  ESCRIPT_DLL_API
  void
  saveVTK(std::string fileName) const;



  ESCRIPT_DLL_API
  Data& operator+=(const Data& right);
  ESCRIPT_DLL_API
  Data& operator+=(const boost::python::object& right);

  ESCRIPT_DLL_API
  Data& operator=(const Data& other);

  ESCRIPT_DLL_API
  Data& operator-=(const Data& right);
  ESCRIPT_DLL_API
  Data& operator-=(const boost::python::object& right);

  ESCRIPT_DLL_API
  Data& operator*=(const Data& right);
  ESCRIPT_DLL_API
  Data& operator*=(const boost::python::object& right);

  ESCRIPT_DLL_API
  Data& operator/=(const Data& right);
  ESCRIPT_DLL_API
  Data& operator/=(const boost::python::object& right);

  ESCRIPT_DLL_API
  Data truedivD(const Data& right);

  ESCRIPT_DLL_API
  Data truedivO(const boost::python::object& right);

  ESCRIPT_DLL_API
  Data rtruedivO(const boost::python::object& left);

  ESCRIPT_DLL_API
  boost::python::object __add__(const boost::python::object& right);
  

  ESCRIPT_DLL_API
  boost::python::object __sub__(const boost::python::object& right);
  
  ESCRIPT_DLL_API
  boost::python::object __rsub__(const boost::python::object& right);  

  ESCRIPT_DLL_API
  boost::python::object __mul__(const boost::python::object& right);
    
  ESCRIPT_DLL_API
  boost::python::object __div__(const boost::python::object& right);
  
  ESCRIPT_DLL_API
  boost::python::object __rdiv__(const boost::python::object& right);    
  
  ESCRIPT_DLL_API
  Data
  matrixInverse() const;

  ESCRIPT_DLL_API
  bool
  probeInterpolation(const FunctionSpace& functionspace) const;

  ESCRIPT_DLL_API
  Data
  getItem(const boost::python::object& key) const;

  ESCRIPT_DLL_API
  void
  setItemD(const boost::python::object& key,
           const Data& value);

  ESCRIPT_DLL_API
  void
  setItemO(const boost::python::object& key,
           const boost::python::object& value);

  // These following public methods should be treated as private.

  template <class UnaryFunction>
  ESCRIPT_DLL_API
  inline
  void
  unaryOp2(UnaryFunction operation);

  ESCRIPT_DLL_API
  Data
  getSlice(const DataTypes::RegionType& region) const;

  ESCRIPT_DLL_API
  void
  setSlice(const Data& value,
           const DataTypes::RegionType& region);

  ESCRIPT_DLL_API
  void
        print(void);

  ESCRIPT_DLL_API
        int
        get_MPIRank(void) const;

  ESCRIPT_DLL_API
        int
        get_MPISize(void) const;

  ESCRIPT_DLL_API
        MPI_Comm
        get_MPIComm(void) const;

  ESCRIPT_DLL_API
        DataAbstract*
        borrowData(void) const;

  ESCRIPT_DLL_API
        DataAbstract_ptr
        borrowDataPtr(void) const;

  ESCRIPT_DLL_API
        DataReady_ptr
        borrowReadyPtr(void) const;



  ESCRIPT_DLL_API
        DataTypes::ValueType::const_reference
        getDataAtOffsetRO(DataTypes::ValueType::size_type i);


  ESCRIPT_DLL_API
        DataTypes::ValueType::reference
        getDataAtOffsetRW(DataTypes::ValueType::size_type i);

 
 protected:

 private:

template <class BinaryOp>
  double 
#ifdef ESYS_MPI
  lazyAlgWorker(double init, MPI_Op mpiop_type);
#else
  lazyAlgWorker(double init);
#endif

  double
  LsupWorker() const;

  double
  supWorker() const;

  double
  infWorker() const;

  boost::python::object
  integrateWorker() const;

  void
  calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;

  void
  calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;

  // For internal use in Data.cpp only!
  // other uses should call the main entry points and allow laziness
  Data
  minval_nonlazy() const;

  // For internal use in Data.cpp only!
  Data
  maxval_nonlazy() const;


  inline
  void
  operandCheck(const Data& right) const
  {
    return m_data->operandCheck(*(right.m_data.get()));
  }

  template <class BinaryFunction>
  inline
  double
  algorithm(BinaryFunction operation,
            double initial_value) const;

  template <class BinaryFunction>
  inline
  Data
  dp_algorithm(BinaryFunction operation,
               double initial_value) const;

  template <class BinaryFunction>
  inline
  void
  binaryOp(const Data& right,
           BinaryFunction operation);

  void
  typeMatchLeft(Data& right) const;

  void
  typeMatchRight(const Data& right);

  void
  initialise(const DataTypes::ValueType& value,
         const DataTypes::ShapeType& shape,
             const FunctionSpace& what,
             bool expanded);

  void
  initialise(const WrappedArray& value,
                 const FunctionSpace& what,
                 bool expanded);

  void
  initialise(const double value,
         const DataTypes::ShapeType& shape,
             const FunctionSpace& what,
             bool expanded);

  //
  // flag to protect the data object against any update
  bool m_protected;
  mutable bool m_shared;
  bool m_lazy;

  //
  // pointer to the actual data object
//   boost::shared_ptr<DataAbstract> m_data;
  DataAbstract_ptr m_data;

// If possible please use getReadyPtr instead.
// But see warning below.
  const DataReady*
  getReady() const;

  DataReady*
  getReady();


// Be wary of using this for local operations since it (temporarily) increases reference count.
// If you are just using this to call a method on DataReady instead of DataAbstract consider using 
// getReady() instead
  DataReady_ptr
  getReadyPtr();

  const_DataReady_ptr
  getReadyPtr() const;


  void updateShareStatus(bool nowshared) const
  {
    m_shared=nowshared;     // m_shared is mutable
  }

  // In the isShared() method below:
  // A problem would occur if m_data (the address pointed to) were being modified 
  // while the call m_data->is_shared is being executed.
  // 
  // Q: So why do I think this code can be thread safe/correct?
  // A: We need to make some assumptions.
  // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
  // 2. We assume that no constructions or assignments which will share previously unshared
  //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
  //
  // This means that the only transition we need to consider, is when a previously shared object is
  // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
  // In those cases the m_shared flag changes to false after m_data has completed changing.
  // For any threads executing before the flag switches they will assume the object is still shared.
  bool isShared() const
  {
    return m_shared;
/*  if (m_shared) return true;
    if (m_data->isShared())         
    {                   
        updateShareStatus(true);
        return true;
    }
    return false;*/
  }

  void forceResolve()
  {
    if (isLazy())
    {
        #ifdef _OPENMP
        if (omp_in_parallel())
        {   // Yes this is throwing an exception out of an omp thread which is forbidden.
        throw DataException("Please do not call forceResolve() in a parallel region.");
        }
        #endif
        resolve();
    }
  }

  void exclusiveWrite()
  {
#ifdef _OPENMP
  if (omp_in_parallel())
  {
// *((int*)0)=17;
    throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
  }
#endif
    forceResolve();
    if (isShared())
    {
        DataAbstract* t=m_data->deepCopy();
        set_m_data(DataAbstract_ptr(t));
    }
  }

  void checkExclusiveWrite()
  {
    if  (isLazy() || isShared())
    {
        throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
    }
  }

  void set_m_data(DataAbstract_ptr p);

  friend class DataAbstract;        // To allow calls to updateShareStatus
  friend class TestDomain;      // so its getX will work quickly
#ifdef IKNOWWHATIMDOING
  friend ESCRIPT_DLL_API Data applyBinaryCFunction(boost::python::object cfunc, boost::python::tuple shape, escript::Data& d, escript::Data& e);
#endif
  friend ESCRIPT_DLL_API Data condEval(escript::Data& mask, escript::Data& trueval, escript::Data& falseval);
  friend ESCRIPT_DLL_API Data randomData(const boost::python::tuple& shape, const FunctionSpace& what, long seed);

};


#ifdef IKNOWWHATIMDOING
ESCRIPT_DLL_API
Data
applyBinaryCFunction(boost::python::object func, boost::python::tuple shape, escript::Data& d, escript::Data& e);
#endif


ESCRIPT_DLL_API
Data
condEval(escript::Data& mask, escript::Data& trueval, escript::Data& falseval);



ESCRIPT_DLL_API
Data randomData(const boost::python::tuple& shape,
       const FunctionSpace& what,
       long seed);


}   // end namespace escript


// No, this is not supposed to be at the top of the file
// DataAbstact needs to be declared first, then DataReady needs to be fully declared
// so that I can dynamic cast between them below.
#include "DataReady.h"
#include "DataLazy.h"

namespace escript
{

inline
const DataReady*
Data::getReady() const
{
   const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
   EsysAssert((dr!=0), "Error - casting to DataReady.");
   return dr;
}

inline
DataReady*
Data::getReady()
{
   DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
   EsysAssert((dr!=0), "Error - casting to DataReady.");
   return dr;
}

// Be wary of using this for local operations since it (temporarily) increases reference count.
// If you are just using this to call a method on DataReady instead of DataAbstract consider using 
// getReady() instead
inline
DataReady_ptr
Data::getReadyPtr()
{
   DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
   EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
   return dr;
}


inline
const_DataReady_ptr
Data::getReadyPtr() const
{
   const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
   EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
   return dr;
}

inline
DataAbstract::ValueType::value_type*
Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
{
   if (isLazy())
   {
    throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
   }
   return getReady()->getSampleDataRW(sampleNo);
}

inline
const DataAbstract::ValueType::value_type*
Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo)
{
   DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
   if (l!=0)
   {
    size_t offset=0;
    const DataTypes::ValueType* res=l->resolveSample(sampleNo,offset);
    return &((*res)[offset]);
   }
   return getReady()->getSampleDataRO(sampleNo);
}



char *Escript_MPI_appendRankToFileName(const char *, int, int);

inline double rpow(double x,double y)
{
    return pow(y,x);
}

ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);

ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);

ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);

ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);

ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);

ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);

ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);

ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);

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);

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);



ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);

ESCRIPT_DLL_API
Data
C_GeneralTensorProduct(Data& arg_0,
                     Data& arg_1,
                     int axis_offset=0,
                     int transpose=0);

inline
Data
Data::truedivD(const Data& right)
{
    return *this / right;
}

inline
Data
Data::truedivO(const boost::python::object& right)
{
    Data tmp(right, getFunctionSpace(), false);
    return truedivD(tmp);
}

inline
Data
Data::rtruedivO(const boost::python::object& left)
{
    Data tmp(left, getFunctionSpace(), false);
    return tmp.truedivD(*this);
}

template <class BinaryFunction>
inline
void
Data::binaryOp(const Data& right,
               BinaryFunction operation)
{
   //
   // if this has a rank of zero promote it to the rank of the RHS
   if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
     throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
   }

   if (isLazy() || right.isLazy())
   {
     throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
   }
   //
   // initially make the temporary a shallow copy
   Data tempRight(right);

   if (getFunctionSpace()!=right.getFunctionSpace()) {
     if (right.probeInterpolation(getFunctionSpace())) {
       //
       // an interpolation is required so create a new Data
       tempRight=Data(right,this->getFunctionSpace());
     } else if (probeInterpolation(right.getFunctionSpace())) {
       //
       // interpolate onto the RHS function space
       Data tempLeft(*this,right.getFunctionSpace());
//        m_data=tempLeft.m_data;
       set_m_data(tempLeft.m_data);
     }
   }
   operandCheck(tempRight);
   //
   // ensure this has the right type for the RHS
   typeMatchRight(tempRight);
   //
   // Need to cast to the concrete types so that the correct binaryOp
   // is called.
   if (isExpanded()) {
     //
     // Expanded data will be done in parallel, the right hand side can be
     // of any data type
     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
     escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
   } else if (isTagged()) {
     //
     // Tagged data is operated on serially, the right hand side can be
     // either DataConstant or DataTagged
     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
     if (right.isTagged()) {
       DataTagged* rightC=dynamic_cast<DataTagged*>(tempRight.m_data.get());
       EsysAssert((rightC!=0), "Programming error - casting to DataTagged.");
       escript::binaryOp(*leftC,*rightC,operation);
     } else {
       DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
       EsysAssert((rightC!=0), "Programming error - casting to DataConstant.");
       escript::binaryOp(*leftC,*rightC,operation);
     }
   } else if (isConstant()) {
     DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
     EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
     escript::binaryOp(*leftC,*rightC,operation);
   }
}

template <class BinaryFunction>
inline
double
Data::algorithm(BinaryFunction operation, double initial_value) const
{
  if (isExpanded()) {
    DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
    EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
    return escript::algorithm(*leftC,operation,initial_value);
  } else if (isTagged()) {
    DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
    EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
    return escript::algorithm(*leftC,operation,initial_value);
  } else if (isConstant()) {
    DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
    EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
    return escript::algorithm(*leftC,operation,initial_value);
  } else if (isEmpty()) {
    throw DataException("Error - Operations not permitted on instances of DataEmpty.");
  } else if (isLazy()) {
    throw DataException("Error - Operations not permitted on instances of DataLazy.");
  } else {
    throw DataException("Error - Data encapsulates an unknown type.");
  }
}

template <class BinaryFunction>
inline
Data
Data::dp_algorithm(BinaryFunction operation, double initial_value) const
{
  if (isEmpty()) {
    throw DataException("Error - Operations not permitted on instances of DataEmpty.");
  } 
  else if (isExpanded()) {
    Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
    DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
    DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
    EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
    EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
    escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
    return result;
  }
  else if (isTagged()) {
    DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
    EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
    DataTypes::ValueType defval(1);
    defval[0]=0;
    DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
    escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
    return Data(resultT);   // note: the Data object now owns the resultT pointer
  } 
  else if (isConstant()) {
    Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
    DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
    DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
    EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
    EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
    escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
    return result;
  } else if (isLazy()) {
    throw DataException("Error - Operations not permitted on instances of DataLazy.");
  } else {
    throw DataException("Error - Data encapsulates an unknown type.");
  }
}

template <typename BinaryFunction>
inline
Data
C_TensorBinaryOperation(Data const &arg_0,
                        Data const &arg_1,
                        BinaryFunction operation)
{
  if (arg_0.isEmpty() || arg_1.isEmpty())
  {
     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
  }
  if (arg_0.isLazy() || arg_1.isLazy())
  {
     throw DataException("Error - Operations not permitted on lazy data.");
  }
  // Interpolate if necessary and find an appropriate function space
  Data arg_0_Z, arg_1_Z;
  if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
    if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
      arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
      arg_1_Z = Data(arg_1);
    }
    else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
      arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
      arg_0_Z =Data(arg_0);
    }
    else {
      throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
    }
  } else {
      arg_0_Z = Data(arg_0);
      arg_1_Z = Data(arg_1);
  }
  // Get rank and shape of inputs
  int rank0 = arg_0_Z.getDataPointRank();
  int rank1 = arg_1_Z.getDataPointRank();
  DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
  DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
  int size0 = arg_0_Z.getDataPointSize();
  int size1 = arg_1_Z.getDataPointSize();
  // Declare output Data object
  Data res;

  if (shape0 == shape1) {
    if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
      double *ptr_2 = &(res.getDataAtOffsetRW(0));

      tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
    }
    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {

      // Prepare the DataConstant input
      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());

      // Borrow DataTagged input from Data object
      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
      res.tag();
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Prepare offset into DataConstant
      int offset_0 = tmp_0->getPointOffset(0,0);
      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));

      // Get the pointers to the actual data
      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));

      // Compute a result for the default
      tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
        tmp_2->addTag(i->first);
        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));

        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
      }

    }
    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_1,dataPointNo_1;
      int numSamples_1 = arg_1_Z.getNumSamples();
      int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
      int offset_0 = tmp_0->getPointOffset(0,0);
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1)); 
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2)); 
          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
      // Borrow DataTagged input from Data object
      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

      // Prepare the DataConstant input
      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
      res.tag();
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Prepare offset into DataConstant
      int offset_1 = tmp_1->getPointOffset(0,0);

      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
      // Get the pointers to the actual data
      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
      // Compute a result for the default
      tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
        tmp_2->addTag(i->first);
        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
      // Borrow DataTagged input from Data object
      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

      // Borrow DataTagged input from Data object
      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
      res.tag();        // DataTagged output
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Get the pointers to the actual data
      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));

      // Compute a result for the default
      tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
      // Merge the tags
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
        tmp_2->addTag(i->first); // use tmp_2 to get correct shape
      }
      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
        tmp_2->addTag(i->first);
      }
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
      for (i=lookup_2.begin();i!=lookup_2.end();i++) {

        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));

        tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      int offset_1 = tmp_1->getPointOffset(0,0);
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);

          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));


          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
      dataPointNo_0=0;
//        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
//       }
      }

    }
    else {
      throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
    }

  } else if (0 == rank0) {
    if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
      double *ptr_2 = &(res.getDataAtOffsetRW(0));
      tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
    }
    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {

      // Prepare the DataConstant input
      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());

      // Borrow DataTagged input from Data object
      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
      res.tag();
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Prepare offset into DataConstant
      int offset_0 = tmp_0->getPointOffset(0,0);
      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));

      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));

      // Compute a result for the default
      tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
        tmp_2->addTag(i->first);
        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
      }

    }
    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {

      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_1;
      int numSamples_1 = arg_1_Z.getNumSamples();
      int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
      int offset_0 = tmp_0->getPointOffset(0,0);
      const double *ptr_src = &(arg_0_Z.getDataAtOffsetRO(offset_0));
      double ptr_0 = ptr_src[0];
      int size = size1*numDataPointsPerSample_1;
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_1) schedule(static)
      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
//        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
          int offset_1 = tmp_1->getPointOffset(sampleNo_1,0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_1,0);
//          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);

//        }
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {

      // Borrow DataTagged input from Data object
      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

      // Prepare the DataConstant input
      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
      res.tag();
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Prepare offset into DataConstant
      int offset_1 = tmp_1->getPointOffset(0,0);
      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));

      // Get the pointers to the actual data
      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));


      // Compute a result for the default
      tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
        tmp_2->addTag(i->first);
        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));

        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {

      // Borrow DataTagged input from Data object
      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

      // Borrow DataTagged input from Data object
      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
      res.tag();        // DataTagged output
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Get the pointers to the actual data
      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));

      // Compute a result for the default
      tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
      // Merge the tags
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
        tmp_2->addTag(i->first); // use tmp_2 to get correct shape
      }
      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
        tmp_2->addTag(i->first);
      }
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));

        tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {

      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      int offset_1 = tmp_1->getPointOffset(0,0);
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
        }
      }


    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {

      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {

      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
        }
      }

    }
    else {
      throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
    }

  } else if (0 == rank1) {
    if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
      double *ptr_2 = &(res.getDataAtOffsetRW(0));
      tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
    }
    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {

      // Prepare the DataConstant input
      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());

      // Borrow DataTagged input from Data object
      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
      res.tag();
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Prepare offset into DataConstant
      int offset_0 = tmp_0->getPointOffset(0,0);
      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));

      //Get the pointers to the actual data
      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));

      // Compute a result for the default
      tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
        tmp_2->addTag(i->first);
        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
      }
    }
    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {

      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_1,dataPointNo_1;
      int numSamples_1 = arg_1_Z.getNumSamples();
      int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
      int offset_0 = tmp_0->getPointOffset(0,0);
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {

      // Borrow DataTagged input from Data object
      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

      // Prepare the DataConstant input
      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
      res.tag();
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Prepare offset into DataConstant
      int offset_1 = tmp_1->getPointOffset(0,0);
      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
      // Get the pointers to the actual data
      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
      // Compute a result for the default
      tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
        tmp_2->addTag(i->first);
        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {

      // Borrow DataTagged input from Data object
      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

      // Borrow DataTagged input from Data object
      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());

      // Prepare a DataTagged output 2
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
      res.tag();        // DataTagged output
      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

      // Get the pointers to the actual data
      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));

      // Compute a result for the default
      tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
      // Merge the tags
      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
        tmp_2->addTag(i->first); // use tmp_2 to get correct shape
      }
      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
        tmp_2->addTag(i->first);
      }
      // Compute a result for each tag
      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
        tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
      }

    }
    else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {

      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      int offset_1 = tmp_1->getPointOffset(0,0);
      const double *ptr_src = &(arg_1_Z.getDataAtOffsetRO(offset_1));
      double ptr_1 = ptr_src[0];
      int size = size0 * numDataPointsPerSample_0;
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
//        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
//          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
//        }
      }


    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {

      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
        }
      }

    }
    else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {

      // After finding a common function space above the two inputs have the same numSamples and num DPPS
      res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

      int sampleNo_0,dataPointNo_0;
      int numSamples_0 = arg_0_Z.getNumSamples();
      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
      res.requireWrite();
      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
        }
      }

    }
    else {
      throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
    }

  } else {
    throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
  }

  return res;
}

template <typename UnaryFunction>
Data
C_TensorUnaryOperation(Data const &arg_0,
                       UnaryFunction operation)
{
  if (arg_0.isEmpty())  // do this before we attempt to interpolate
  {
     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
  }
  if (arg_0.isLazy())
  {
     throw DataException("Error - Operations not permitted on lazy data.");
  }
  // Interpolate if necessary and find an appropriate function space
  Data arg_0_Z = Data(arg_0);

  // Get rank and shape of inputs
  const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
  int size0 = arg_0_Z.getDataPointSize();

  // Declare output Data object
  Data res;

  if (arg_0_Z.isConstant()) {
    res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
    const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
    double *ptr_2 = &(res.getDataAtOffsetRW(0));
    tensor_unary_operation(size0, ptr_0, ptr_2, operation);
  }
  else if (arg_0_Z.isTagged()) {

    // Borrow DataTagged input from Data object
    DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());

    // Prepare a DataTagged output 2
    res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
    res.tag();
    DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());

    // Get the pointers to the actual data
    const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
    double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
    // Compute a result for the default
    tensor_unary_operation(size0, ptr_0, ptr_2, operation);
    // Compute a result for each tag
    const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
    DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
    for (i=lookup_0.begin();i!=lookup_0.end();i++) {
      tmp_2->addTag(i->first);
      const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
      double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
      tensor_unary_operation(size0, ptr_0, ptr_2, operation);
    }

  }
  else if (arg_0_Z.isExpanded()) {

    res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
    DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
    DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());

    int sampleNo_0,dataPointNo_0;
    int numSamples_0 = arg_0_Z.getNumSamples();
    int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
    #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
    for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
    dataPointNo_0=0;
//      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
        int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
        int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
        double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
        tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
//      }
    }
  }
  else {
    throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
  }

  return res;
}

}
#endif