SharedArray.hpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File SharedArray.hpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Scientific Computing and Imaging Institute,
// University of Utah (USA) and Department of Aeronautics, Imperial
// College London (UK).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description:
//
///////////////////////////////////////////////////////////////////////////////
 
#ifndef NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHARED_ARRAY_HPP
#define NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHARED_ARRAY_HPP
 
#include <LibUtilities/BasicUtils/ArrayPolicies.hpp>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/LinearAlgebra/NekVectorFwd.hpp>
#include <LibUtilities/LinearAlgebra/NekMatrixFwd.hpp>
#include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
 
#include <boost/core/ignore_unused.hpp>
#include <boost/multi_array.hpp>
 
namespace Nektar
{
    class LinearSystem;
 
    // Forward declaration for a ConstArray constructor.
    template<typename Dim, typename DataType>
    class Array;
 
    /// \brief 1D Array of constant elements with garbage collection and bounds checking.
    template<typename DataType>
    class Array<OneD, const DataType>
    {
#ifdef WITH_PYTHON
        struct PythonInfo {
            void *m_pyObject; // Underlying PyObject pointer
            void (*m_callback)(void *); // Callback
        };
#endif
        public:
            typedef DataType* ArrayType;
            typedef const DataType& const_reference;
            typedef DataType& reference;
 
            typedef const DataType* const_iterator;
            typedef DataType* iterator;
 
            typedef DataType element;
            typedef size_t size_type;
 
 
        public:
            /// \brief Creates an empty array.
            Array() :
#ifdef WITH_PYTHON
                m_pythonInfo(nullptr),
#endif
                m_size( 0 ),
                m_capacity( 0 ),
                m_data( nullptr ),
                m_count( nullptr ),
                m_offset( 0 )                
            {
                CreateStorage(m_capacity);
            }
 
            /// \brief Creates an array of size dim1Size.
            ///
            /// If DataType is a fundamental type (double, int, etc.), then the allocated array is
            /// uninitialized.  If it is any other type, each element is initialized with DataType's default
            /// constructor.
            explicit Array(size_t dim1Size) :
#ifdef WITH_PYTHON
                m_pythonInfo(nullptr),
#endif
                m_size( dim1Size ),
                m_capacity( dim1Size ),
                m_data( nullptr ),
                m_count( nullptr ),
                m_offset( 0 )                
            {
                CreateStorage(m_capacity);
                ArrayInitializationPolicy<DataType>::Initialize( m_data, m_capacity );
            }
 
            /// \brief Creates a 1D array with each element
            /// initialized to an initial value.
            /// \param dim1Size The array's size.
            /// \param initValue Each element's initial value.
            ///
            /// If DataType is a fundamental type (double, int, etc.),
            /// then the initial value is copied directly into each
            /// element.  Otherwise, the DataType's copy constructor
            /// is used to initialize each element.
            Array(size_t dim1Size, const DataType& initValue) :
#ifdef WITH_PYTHON
                m_pythonInfo(nullptr),
#endif
                m_size( dim1Size ),
                m_capacity( dim1Size ),
                m_data( nullptr ),
                m_count( nullptr ),
                m_offset( 0 )                
            {
                CreateStorage(m_capacity);
                ArrayInitializationPolicy<DataType>::Initialize( m_data, m_capacity, initValue );
            }
 
            /// \brief Creates a 1D array a copies data into it.
            /// \param dim1Size the array's size.
            /// \param data The data to copy.
            ///
            /// If DataType is a fundamental type (double, int, etc.), then data is copied
            /// directly into the underlying storage.  Otherwise, the DataType's copy constructor
            /// is used to copy each element.
            Array(size_t dim1Size, const DataType* data) :
#ifdef WITH_PYTHON
                m_pythonInfo(nullptr),
#endif
                m_size( dim1Size ),
                m_capacity( dim1Size ),
                m_data( nullptr ),
                m_count( nullptr ),
                m_offset( 0 )                
            {
                CreateStorage(m_capacity);
                ArrayInitializationPolicy<DataType>::Initialize( m_data, m_capacity, data );
            }
 
            /// \brief Creates a 1D array that references rhs.
            /// \param dim1Size The size of the array.  This is useful
            ///                 when you want this array to reference
            ///                 a subset of the elements in rhs.
            /// \param rhs      Array to reference.
            /// This constructor creates an array that references rhs.
            /// Any changes to rhs will be reflected in this array.
            /// The memory for the array will only be deallocated when
            /// both rhs and this array have gone out of scope.
            Array(size_t dim1Size, const Array<OneD, const DataType>& rhs) :
#ifdef WITH_PYTHON
                m_pythonInfo(rhs.m_pythonInfo),
#endif
                m_size(dim1Size),
                m_capacity(rhs.m_capacity),
                m_data(rhs.m_data),
                m_count(rhs.m_count),
                m_offset(rhs.m_offset)                
            {
                *m_count += 1;
                ASSERTL0(m_size <= rhs.num_elements(), "Requested size is larger than input array size.");
            }
 
#ifdef WITH_PYTHON
            /// \brief Creates a 1D array a copies data into it.
            /// \param dim1Size the array's size.
            /// \param data The data to reference.
            /// \param memory_pointer Pointer to the memory address of the array
            /// \param python_decrement Pointer to decrementer
            Array(size_t dim1Size, DataType* data, void* memory_pointer, void (*python_decrement)(void *)) :
                m_size( dim1Size ),
                m_capacity( dim1Size ),
                m_data( data ),
                m_count( nullptr ),
                m_offset( 0 )                               
            {
                m_count = new size_t(); 
                *m_count = 1;
 
                m_pythonInfo = new PythonInfo *();
                *m_pythonInfo = new PythonInfo();
                (*m_pythonInfo)->m_callback = python_decrement;
                (*m_pythonInfo)->m_pyObject = memory_pointer;
            }
#endif
 
            /// \brief Creates a reference to rhs.
            Array(const Array<OneD, const DataType>& rhs) :
#ifdef WITH_PYTHON
                m_pythonInfo(rhs.m_pythonInfo),
#endif
                m_size(rhs.m_size),
                m_capacity(rhs.m_capacity),
                m_data(rhs.m_data),
                m_count(rhs.m_count),
                m_offset(rhs.m_offset)                
            {
                *m_count += 1;
            }
 
            ~Array()
            {
                if( m_count == nullptr )
                {
                    return;
                }
 
                *m_count -= 1;
                if( *m_count == 0 )
                {
#ifdef WITH_PYTHON
                    if (*m_pythonInfo == nullptr)
                    {
                        ArrayDestructionPolicy<DataType>::Destroy( m_data, m_capacity );
                        MemoryManager<DataType>::RawDeallocate( m_data, m_capacity );
                    }
                    else
                    {
                        (*m_pythonInfo)->m_callback((*m_pythonInfo)->m_pyObject);
                        delete *m_pythonInfo;
                    }
 
                    delete m_pythonInfo;
 
#else
 
                    ArrayDestructionPolicy<DataType>::Destroy( m_data, m_capacity );
                    MemoryManager<DataType>::RawDeallocate( m_data, m_capacity );
 
#endif
 
                    delete m_count; // Clean up the memory used for the reference count.
                }
            }
 
            /// \brief Creates a reference to rhs.
            Array<OneD, const DataType>& operator=(const Array<OneD, const DataType>& rhs) 
            {
                *m_count -= 1;
                if( *m_count == 0 )
                {
#ifdef WITH_PYTHON
                    if (*m_pythonInfo == nullptr)
                    {
                        ArrayDestructionPolicy<DataType>::Destroy( m_data, m_capacity );
                        MemoryManager<DataType>::RawDeallocate( m_data, m_capacity );
                    }
                    else if ((*rhs.m_pythonInfo) != nullptr && (*m_pythonInfo)->m_pyObject != (*rhs.m_pythonInfo)->m_pyObject)
                    {
                        (*m_pythonInfo)->m_callback((*m_pythonInfo)->m_pyObject);
                        delete *m_pythonInfo;
                    }
 
                    delete m_pythonInfo;
#else
 
                    ArrayDestructionPolicy<DataType>::Destroy( m_data, m_capacity );
                    MemoryManager<DataType>::RawDeallocate( m_data, m_capacity );
#endif
                    delete m_count; // Clean up the memory used for the reference count.
                }
 
                m_data = rhs.m_data;
                m_capacity = rhs.m_capacity;
                m_count = rhs.m_count;
                *m_count += 1;
                m_offset = rhs.m_offset;
                m_size = rhs.m_size;
#ifdef WITH_PYTHON
                m_pythonInfo = rhs.m_pythonInfo;
#endif
                return *this;
            }
 
            const_iterator begin() const { return m_data + m_offset; }
            const_iterator end() const { return m_data + m_offset + m_size; }
 
            const_reference operator[](size_t i) const
            {
                ASSERTL1(i < m_size,
                         std::string("Element ") + std::to_string(i) +
                         std::string(" requested in an array of size ") +
                         std::to_string(m_size));
                return *( m_data + i + m_offset );
            }
 
            /// \brief Returns a c-style pointer to the underlying array.
            const element* get() const { return m_data + m_offset; }
 
            /// \brief Returns a c-style pointer to the underlying array.
            const element* data() const { return m_data + m_offset; }
 
            /// \brief Returns 1.
            size_t num_dimensions() const { return 1; }
 
            /// \brief Returns the array's size.
            size_t num_elements() const { return m_size; }
 
            size_t capacity() const { return m_capacity; }
 
            /// \brief Returns the array's offset.
            size_t GetOffset() const { return m_offset; }
 
            /// \brief Returns the array's reference counter.
            size_t GetCount() const { return *m_count; }
 
            /// \brief Returns true is this array and rhs overlap.
            bool Overlaps(const Array<OneD, const DataType>& rhs) const
            {
                const element* start = get();
                const element* end = start + m_size;
 
                const element* rhs_start = rhs.get();
                const element* rhs_end = rhs_start + rhs.num_elements();
 
                return (rhs_start >= start && rhs_start <= end) ||
                       (rhs_end >= start && rhs_end <= end);
            }
 
#ifdef WITH_PYTHON
            bool IsPythonArray()
            {
                return *m_pythonInfo != nullptr;
            }
 
            void ToPythonArray(void* memory_pointer, void (*python_decrement)(void *))
            {
                *m_pythonInfo = new PythonInfo();
                (*m_pythonInfo)->m_callback = python_decrement;
                (*m_pythonInfo)->m_pyObject = memory_pointer;
            }
#endif
 
            template<typename T1, typename T2>
            friend bool operator==(const Array<OneD, T1>&, const Array<OneD, T2>&);
            LIB_UTILITIES_EXPORT friend bool operator==(const Array<OneD, NekDouble>&, const Array<OneD, NekDouble>&);
            LIB_UTILITIES_EXPORT friend bool IsEqual(const Array<OneD, const NekDouble>&,const Array<OneD, const NekDouble>&,NekDouble);
 
            /// \brief Creates an array with a specified offset.
            ///
            /// The return value will reference the same array as lhs,
            /// but with an offset.
            ///
            /// For example, in the following:
            /// \code
            /// Array<OneD, const double> result = anArray + 10;
            /// \endcode
            /// result[0] == anArray[10];
            template<typename T>
            friend Array<OneD, T> operator+(const Array<OneD, T>& lhs, size_t offset);
 
            template<typename T>
            friend Array<OneD, T> operator+(unsigned int offset, const Array<OneD, T>& rhs);
 
        protected:
 
#ifdef WITH_PYTHON
            PythonInfo **m_pythonInfo;
#endif
 
            size_t m_size;
            size_t m_capacity;
            DataType* m_data;
 
            // m_count points to an integer used as a reference count to this array's data (m_data).
            // Previously, the reference count was stored in the first 4 bytes of the m_data array.
            size_t* m_count;
 
            size_t m_offset;
 
        private:
        //            struct DestroyArray
        //            {
        //                DestroyArray(unsigned int elements) :
        //                    m_elements(elements) {}
        //
        //                void operator()(DataType* p)
        //                {
        //                    ArrayDestructionPolicy<DataType>::Destroy(p, m_elements);
        //                    MemoryManager<DataType>::RawDeallocate(p, m_elements);
        //                }
        //                unsigned int m_elements;
        //            };
        //
        void
            CreateStorage( size_t size )
            {
                DataType* storage = MemoryManager<DataType>::RawAllocate( size );
                m_data = storage;
 
                // Allocate an integer to hold the reference count.  Note 1, all arrays that share this array's
                // data (ie, point to m_data) will also share the m_count data.  Note 2, previously m_count
                // pointed to "(unsigned int*)storage".
                m_count = new size_t();
                *m_count = 1;
#ifdef WITH_PYTHON
                m_pythonInfo = new PythonInfo*();
                *m_pythonInfo = nullptr;
#endif
            }
 
            template<typename T>
            static Array<OneD, T> CreateWithOffset(const Array<OneD, T>& rhs, size_t offset)
            {
                Array<OneD, T> result(rhs);
                result.m_offset += offset;
                result.m_size = rhs.m_size - offset;
                return result;
            }
 
    };
 
 
    /// \brief 2D array with garbage collection and bounds checking.
    template<typename DataType>
    class Array<TwoD, const DataType>
    {
        public:
            typedef boost::multi_array_ref<DataType, 2> ArrayType;
            typedef typename ArrayType::const_reference const_reference;
            typedef typename ArrayType::reference reference;
 
            typedef typename ArrayType::index index;
            typedef typename ArrayType::const_iterator const_iterator;
            typedef typename ArrayType::iterator iterator;
 
            typedef typename ArrayType::element element;
            typedef typename ArrayType::size_type size_type;
 
 
 
        public:
            Array() :
                m_data(CreateStorage<DataType>(0, 0))
            {
            }
 
            /// \brief Constructs a 2 dimensional array.  The elements of the array are not initialized.
            Array(size_t dim1Size, size_t dim2Size) :
                m_data(CreateStorage<DataType>(dim1Size, dim2Size))
            {
                ArrayInitializationPolicy<DataType>::Initialize(m_data->data(), m_data->num_elements());
            }
 
            Array(size_t dim1Size, size_t dim2Size, const DataType& initValue) :
                m_data(CreateStorage<DataType>(dim1Size, dim2Size))
            {
                ArrayInitializationPolicy<DataType>::Initialize(m_data->data(), m_data->num_elements(), initValue);
            }
 
            Array(size_t dim1Size, size_t dim2Size, const DataType* data) :
                m_data(CreateStorage<DataType>(dim1Size, dim2Size))
            {
                ArrayInitializationPolicy<DataType>::Initialize(m_data->data(), m_data->num_elements(), data);
            }
 
            Array(const Array<TwoD, const DataType>& rhs) :
                m_data(rhs.m_data)
            {
            }
 
            Array<TwoD, const DataType>& operator=(const Array<TwoD, const DataType>& rhs)
            {
                m_data = rhs.m_data;
                return *this;
            }
 
            template<typename T>
            friend bool operator==(const Array<TwoD, T>&, const Array<TwoD, T>&);
            LIB_UTILITIES_EXPORT friend bool operator==(const Array<TwoD, NekDouble>&, const Array<TwoD, NekDouble>&);
            LIB_UTILITIES_EXPORT friend bool IsEqual(const Array<TwoD, const NekDouble>&,const Array<TwoD, const NekDouble>&,NekDouble);
 
            const_iterator begin() const { return m_data->begin(); }
            const_iterator end() const { return m_data->end(); }
            const_reference operator[](index i) const { return (*m_data)[i]; }
            const element* get() const { return m_data->data(); }
            const element* data() const { return m_data->data(); }
            size_t num_dimensions() const { return m_data->num_dimensions(); }
            const size_type* shape() const { return m_data->shape(); }
            size_t num_elements() const { return m_data->num_elements(); }
 
            size_t GetRows() const { return m_data->shape()[0]; }
            size_t GetColumns() const { return m_data->shape()[1]; }
 
        protected:
            std::shared_ptr<ArrayType> m_data;
 
        private:
 
    };
 
    enum AllowWrappingOfConstArrays
    {
        eVECTOR_WRAPPER
    };
 
    /// \brief 1D Array
    ///
    /// \ref pageNekArrays
    ///
    /// Misc notes.
    ///
    /// Through out the 1D Array class you will see things like "using BaseType::begin" and
    /// "using BaseType::end".  This is necessary to bring the methods from the ConstArray
    /// into scope in Array class.  Typically this is not necessary, but since we have
    /// method names which match those in the base class, the base class names are hidden.
    /// Therefore, we have to explicitly bring them into scope to use them.
    template<typename DataType>
    class Array<OneD, DataType> : public Array<OneD, const DataType>
    {
        public:
            typedef Array<OneD, const DataType> BaseType;
            typedef typename BaseType::iterator iterator;
            typedef typename BaseType::reference reference;
            typedef typename BaseType::size_type size_type;
            typedef typename BaseType::element element;
 
        public:
            Array() :
                BaseType()
            {
            }
 
            explicit Array(size_t dim1Size) :
                BaseType(dim1Size)
            {
            }
 
            Array(size_t dim1Size, const DataType& initValue) :
                BaseType(dim1Size, initValue)
            {
            }
 
            Array(size_t dim1Size, const DataType* data) :
                BaseType(dim1Size, data)
            {
            }
 
            Array(size_t dim1Size, const Array<OneD, DataType>& rhs) :
                BaseType(dim1Size, rhs)
            {
            }
 
            Array(size_t dim1Size, const Array<OneD, const DataType>& rhs) :
                BaseType(dim1Size, rhs.data())
            {
            }
 
            Array(const Array<OneD, DataType>& rhs) :
                BaseType(rhs)
            {
            }
 
            Array(const Array<OneD, const DataType>& rhs) :
                BaseType(rhs.num_elements(), rhs.data())
            {
            }
 
#ifdef WITH_PYTHON
            Array(size_t dim1Size, DataType* data, void* memory_pointer, void (*python_decrement)(void *)) :
                BaseType(dim1Size, data, memory_pointer, python_decrement)
            {
            }
 
            void ToPythonArray(void* memory_pointer, void (*python_decrement)(void *))
            {
                BaseType::ToPythonArray(memory_pointer, python_decrement);
            }
#endif
 
            Array<OneD, DataType>& operator=(const Array<OneD, DataType>& rhs)
            {
                BaseType::operator=(rhs);
                return *this;
            }
 
            static Array<OneD, DataType> CreateWithOffset(const Array<OneD, DataType>& rhs, unsigned int offset)
            {
                Array<OneD, DataType> result(rhs);
                result.m_offset += offset;
                result.m_size = rhs.m_size - offset;
                return result;
            }
 
            typename Array<OneD, const DataType>::const_iterator
                     begin() const { return Array<OneD, const DataType>::begin(); }
            iterator begin() { return this->m_data + this->m_offset; }
 
            using BaseType::end;
            iterator end() { return this->m_data + this->m_offset + this->m_size; }
 
            using BaseType::operator[];
            reference operator[](size_t i)
            {
                ASSERTL1(static_cast<size_type>(i) < this->num_elements(),
                         std::string("Element ") + std::to_string(i) +
                         std::string(" requested in an array of size ") +
                         std::to_string(this->num_elements()));
                return (get())[i];
            }
 
 
            using BaseType::get;
            element* get() { return this->m_data + this->m_offset; }
 
            using BaseType::data;
            element* data() { return this->m_data + this->m_offset; }
 
            template<typename T1>
            friend class NekVector;
 
            template<typename T1, typename T3>
            friend class NekMatrix;
 
            friend class LinearSystem;
 
        protected:
            Array(const Array<OneD, const DataType>& rhs, AllowWrappingOfConstArrays a) :
                BaseType(rhs)
            {
                boost::ignore_unused(a);
            }
 
            void ChangeSize(size_t newSize)
            {
                ASSERTL1(newSize <= this->m_capacity, "Can't change an array size to something larger than its capacity.");
                this->m_size = newSize;
            }
 
        private:
 
    };
 
    /// \brief A 2D array.
    template<typename DataType>
    class Array<TwoD, DataType> : public Array<TwoD, const DataType>
    {
        public:
            typedef Array<TwoD, const DataType> BaseType;
            typedef typename BaseType::iterator iterator;
            typedef typename BaseType::reference reference;
            typedef typename BaseType::index index;
            typedef typename BaseType::size_type size_type;
            typedef typename BaseType::element element;
 
        public:
            Array() :
                BaseType()
            {
            }
 
            Array(size_t dim1Size, size_t dim2Size) :
                BaseType(dim1Size, dim2Size)
            {
            }
 
            Array(size_t dim1Size, size_t dim2Size, const DataType& initValue) :
                BaseType(dim1Size, dim2Size, initValue)
            {
            }
 
            Array(size_t dim1Size, size_t dim2Size, const DataType* data) :
                BaseType(dim1Size, dim2Size, data)
            {
            }
 
            Array(const Array<TwoD, DataType>& rhs) :
                BaseType(rhs)
            {
            }
 
            Array<TwoD, DataType>& operator=(const Array<TwoD, DataType>& rhs)
            {
                BaseType::operator=(rhs);
                return *this;
            }
 
            using BaseType::begin;
            iterator begin() { return this->m_data->begin(); }
 
            using BaseType::end;
            iterator end() { return this->m_data->end(); }
 
            using BaseType::operator[];
            reference operator[](index i) { return (*this->m_data)[i]; }
 
            using BaseType::get;
            element* get() { return this->m_data->data(); }
 
            using BaseType::data;
            element* data() { return this->m_data->data(); }
 
        private:
 
    };
 
    LIB_UTILITIES_EXPORT bool IsEqual(const Array<OneD, const NekDouble>& lhs,
                 const Array<OneD, const NekDouble>& rhs,
                 NekDouble tol = NekConstants::kNekZeroTol);
    LIB_UTILITIES_EXPORT bool operator==(const Array<OneD, NekDouble>& lhs, const Array<OneD, NekDouble>& rhs);
 
    template<typename T1, typename T2>
    bool operator==(const Array<OneD, T1>& lhs, const Array<OneD, T2>& rhs)
    {
        if( lhs.num_elements() != rhs.num_elements() )
        {
            return false;
        }
 
        if( lhs.data() == rhs.data() )
        {
            return true;
        }
 
        for(unsigned int i = 0; i < lhs.num_elements(); ++i)
        {
            if( lhs[i] != rhs[i] )
            {
                return false;
            }
        }
 
        return true;
    }
 
    template<typename T1, typename T2>
    bool operator!=(const Array<OneD, T1>& lhs, const Array<OneD, T2>& rhs)
    {
        return !(lhs == rhs);
    }
 
    template<typename DataType>
    Array<OneD, DataType> operator+(const Array<OneD, DataType>& lhs, size_t offset)
    {
        return Array<OneD, const DataType>::CreateWithOffset(lhs, offset);
    }
 
    template<typename DataType>
    Array<OneD, DataType> operator+(size_t offset, const Array<OneD, DataType>& rhs)
    {
        return Array<OneD, const DataType>::CreateWithOffset(rhs, offset);
    }
 
//    template<typename DataType>
//    Array<OneD, DataType> operator+(const Array<OneD, DataType>& lhs, size_t offset)
//    {
//        return Array<OneD, DataType>::CreateWithOffset(lhs, offset);
//    }
 
    template<typename ConstDataType, typename DataType>
    void CopyArray(const Array<OneD, ConstDataType>& source, Array<OneD, DataType>& dest)
    {
        if( dest.num_elements() != source.num_elements() )
        {
            dest = Array<OneD, DataType>(source.num_elements());
        }
 
        std::copy(source.data(), source.data() + source.num_elements(), dest.data());
    }
 
    template<typename ConstDataType, typename DataType>
    void CopyArrayN(const Array<OneD, ConstDataType>& source, Array<OneD, DataType>& dest, unsigned int n)
    {
        if( dest.num_elements() != n )
        {
            dest = Array<OneD, DataType>(n);
        }
 
        std::copy(source.data(), source.data() + n, dest.data());
    }
 
    static Array<OneD, int> NullInt1DArray;
    static Array<OneD, NekDouble> NullNekDouble1DArray;
    static Array<OneD, Array<OneD, NekDouble> > NullNekDoubleArrayofArray;
 
    LIB_UTILITIES_EXPORT bool IsEqual(const Array<TwoD, const NekDouble>& lhs,
                 const Array<TwoD, const NekDouble>& rhs,
                 NekDouble tol = NekConstants::kNekZeroTol);
    LIB_UTILITIES_EXPORT bool operator==(const Array<TwoD, NekDouble>& lhs, const Array<TwoD, NekDouble>& rhs) ;
 
    template<typename DataType>
    bool operator==(const Array<TwoD, DataType>& lhs, const Array<TwoD, DataType>& rhs)
    {
        return *lhs.m_data == *rhs.m_data;
    }
 
    template<typename DataType>
    bool operator!=(const Array<TwoD, DataType>& lhs, const Array<TwoD, DataType>& rhs)
    {
        return !(lhs == rhs);
    }
 
    namespace LibUtilities
    {
        static std::vector<NekDouble> NullNekDoubleVector;
        static std::vector<unsigned int> NullUnsignedIntVector;
        static std::vector<std::vector<NekDouble> > NullVectorNekDoubleVector
            = { NullNekDoubleVector };
    }
}
 
#endif //NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHARED_ARRAY_HPP