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