diff --git a/src/TNL/Containers/Algorithms/ArrayOperations.h b/src/TNL/Containers/Algorithms/ArrayOperations.h
index a6acc816edc4e20a73949c49df42164e365475b9..c7a8fe09424066110ac8a1a5cf652fb48daa4769 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperations.h
+++ b/src/TNL/Containers/Algorithms/ArrayOperations.h
@@ -28,13 +28,6 @@ class ArrayOperations< Devices::Host >
 {
    public:
 
-      template< typename Element, typename Index >
-      static void allocateMemory( Element*& data,
-                                  const Index size );
-
-      template< typename Element >
-      static void freeMemory( Element* data );
-
       template< typename Element >
       static void setMemoryElement( Element* data,
                                     const Element& value );
@@ -84,13 +77,6 @@ class ArrayOperations< Devices::Cuda >
 {
    public:
 
-      template< typename Element, typename Index >
-      static void allocateMemory( Element*& data,
-                                  const Index size );
-
-      template< typename Element >
-      static void freeMemory( Element* data );
-
       template< typename Element >
       static void setMemoryElement( Element* data,
                                     const Element& value );
@@ -181,13 +167,6 @@ class ArrayOperations< Devices::MIC >
 {
    public:
 
-      template< typename Element, typename Index >
-      static void allocateMemory( Element*& data,
-                                  const Index size );
-
-      template< typename Element >
-      static void freeMemory( Element* data );
-
       template< typename Element >
       static void setMemoryElement( Element* data,
                                     const Element& value );
diff --git a/src/TNL/Containers/Algorithms/ArrayOperationsCuda.hpp b/src/TNL/Containers/Algorithms/ArrayOperationsCuda.hpp
index 3e696672c9ef1d4f817db36253dc3e28ecb47146..fe5637c4d9bd2cd7e93e6429cd14b8f664dcc18c 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperationsCuda.hpp
+++ b/src/TNL/Containers/Algorithms/ArrayOperationsCuda.hpp
@@ -16,7 +16,6 @@
 #include <TNL/Math.h>
 #include <TNL/ParallelFor.h>
 #include <TNL/Exceptions/CudaSupportMissing.h>
-#include <TNL/Exceptions/CudaBadAlloc.h>
 #include <TNL/Containers/Algorithms/ArrayOperations.h>
 #include <TNL/Containers/Algorithms/Reduction.h>
 #include <TNL/Containers/Algorithms/ReductionOperations.h>
@@ -25,41 +24,6 @@ namespace TNL {
 namespace Containers {
 namespace Algorithms {
 
-template< typename Element, typename Index >
-void
-ArrayOperations< Devices::Cuda >::
-allocateMemory( Element*& data,
-                const Index size )
-{
-#ifdef HAVE_CUDA
-   TNL_CHECK_CUDA_DEVICE;
-   if( cudaMalloc( ( void** ) &data,
-                   ( std::size_t ) size * sizeof( Element ) ) != cudaSuccess )
-   {
-      data = 0;
-      throw Exceptions::CudaBadAlloc();
-   }
-   TNL_CHECK_CUDA_DEVICE;
-#else
-   throw Exceptions::CudaSupportMissing();
-#endif
-}
-
-template< typename Element >
-void
-ArrayOperations< Devices::Cuda >::
-freeMemory( Element* data )
-{
-   TNL_ASSERT_TRUE( data, "Attempted to free a nullptr." );
-#ifdef HAVE_CUDA
-   TNL_CHECK_CUDA_DEVICE;
-   cudaFree( data );
-   TNL_CHECK_CUDA_DEVICE;
-#else
-   throw Exceptions::CudaSupportMissing();
-#endif
-}
-
 template< typename Element >
 void
 ArrayOperations< Devices::Cuda >::
diff --git a/src/TNL/Containers/Algorithms/ArrayOperationsHost.hpp b/src/TNL/Containers/Algorithms/ArrayOperationsHost.hpp
index 390b9120754ee6d16f9f551cbab691f45fc102cc..6f4751d4a31f05914fc5b4afcb48e085d20b060d 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperationsHost.hpp
+++ b/src/TNL/Containers/Algorithms/ArrayOperationsHost.hpp
@@ -22,28 +22,6 @@ namespace TNL {
 namespace Containers {
 namespace Algorithms {
 
-template< typename Element, typename Index >
-void
-ArrayOperations< Devices::Host >::
-allocateMemory( Element*& data,
-                const Index size )
-{
-   data = new Element[ size ];
-   // According to the standard, new either throws, or returns non-nullptr.
-   // Some (old) compilers don't comply:
-   // https://stackoverflow.com/questions/550451/will-new-return-null-in-any-case
-   TNL_ASSERT_TRUE( data, "Operator 'new' returned a nullptr. This should never happen - there is "
-                          "either a bug or the compiler does not comply to the standard." );
-}
-
-template< typename Element >
-void
-ArrayOperations< Devices::Host >::
-freeMemory( Element* data )
-{
-   delete[] data;
-}
-
 template< typename Element >
 void
 ArrayOperations< Devices::Host >::
diff --git a/src/TNL/Containers/Algorithms/ArrayOperationsMIC.hpp b/src/TNL/Containers/Algorithms/ArrayOperationsMIC.hpp
index b0f09fb026fcdb2f0124cba6657b24c269c8c96a..fc0e7f56b03d7a33c2becce56149ae8e38756e4f 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperationsMIC.hpp
+++ b/src/TNL/Containers/Algorithms/ArrayOperationsMIC.hpp
@@ -16,7 +16,6 @@
 
 #include <TNL/Math.h>
 #include <TNL/Exceptions/MICSupportMissing.h>
-#include <TNL/Exceptions/MICBadAlloc.h>
 #include <TNL/Containers/Algorithms/ArrayOperations.h>
 #include <TNL/Containers/Algorithms/Reduction.h>
 #include <TNL/Containers/Algorithms/ReductionOperations.h>
@@ -28,34 +27,6 @@ namespace Algorithms {
 
 static constexpr std::size_t MIC_STACK_VAR_LIM = 5*1024*1024;
 
-template< typename Element, typename Index >
-void
-ArrayOperations< Devices::MIC >::
-allocateMemory( Element*& data,
-                const Index size )
-{
-#ifdef HAVE_MIC
-   data = (Element*) Devices::MIC::AllocMIC( size * sizeof(Element) );
-   if( ! data )
-      throw Exceptions::MICBadAlloc();
-#else
-   throw Exceptions::MICSupportMissing();
-#endif
-}
-
-template< typename Element >
-void
-ArrayOperations< Devices::MIC >::
-freeMemory( Element* data )
-{
-   TNL_ASSERT( data, );
-#ifdef HAVE_MIC
-   Devices::MIC::FreeMIC( data );
-#else
-   throw Exceptions::MICSupportMissing();
-#endif
-}
-
 template< typename Element >
 void
 ArrayOperations< Devices::MIC >::
diff --git a/src/TNL/Containers/Array.h b/src/TNL/Containers/Array.h
index 06cae8787ffeec6097070a69f013de29c22eea32..9daa3e2d27afb09852234f987f6cfe1221185fff 100644
--- a/src/TNL/Containers/Array.h
+++ b/src/TNL/Containers/Array.h
@@ -15,6 +15,7 @@
 
 #include <TNL/File.h>
 #include <TNL/TypeTraits.h>
+#include <TNL/Allocators/Default.h>
 #include <TNL/Containers/ArrayView.h>
 
 namespace TNL {
@@ -34,6 +35,8 @@ template< int, typename > class StaticArray;
  *                the compile-time checks of correct pointers manipulation. It
  *                can be either \ref Devices::Host or \ref Devices::Cuda.
  * \tparam Index  The indexing type.
+ * \tparam Allocator The type of the allocator used for the allocation and
+ *                   deallocation of memory used by the array.
  *
  * Memory management handled by constructors and destructors according to the
  * [RAII](https://en.wikipedia.org/wiki/RAII) principle and by methods
@@ -64,7 +67,8 @@ template< int, typename > class StaticArray;
  */
 template< typename Value,
           typename Device = Devices::Host,
-          typename Index = int >
+          typename Index = int,
+          typename Allocator = typename Allocators::Default< Device >::template Allocator< Value > >
 class Array
 {
    public:
@@ -72,6 +76,7 @@ class Array
       using ValueType = Value;
       using DeviceType = Device;
       using IndexType = Index;
+      using AllocatorType = Allocator;
       using HostType = Containers::Array< Value, Devices::Host, Index >;
       using CudaType = Containers::Array< Value, Devices::Cuda, Index >;
       using ViewType = ArrayView< Value, Device, Index >;
@@ -80,14 +85,22 @@ class Array
       /**
        * \brief Constructs an empty array with zero size.
        */
-      Array();
+      Array() = default;
+
+      /**
+       * \brief Constructs an empty array and sets the provided allocator.
+       *
+       * \param allocator The allocator to be associated with this array.
+       */
+      explicit Array( const AllocatorType& allocator );
 
       /**
        * \brief Constructs an array with given size.
        *
        * \param size The number of array elements to be allocated.
+       * \param allocator The allocator to be associated with this array.
        */
-      Array( const IndexType& size );
+      explicit Array( const IndexType& size, const AllocatorType& allocator = AllocatorType() );
 
       /**
        * \brief Constructs an array with given size and copies data from given
@@ -95,9 +108,11 @@ class Array
        *
        * \param data The pointer to the data to be copied to the array.
        * \param size The number of array elements to be copied to the array.
+       * \param allocator The allocator to be associated with this array.
        */
       Array( Value* data,
-             const IndexType& size );
+             const IndexType& size,
+             const AllocatorType& allocator = AllocatorType() );
 
       /**
        * \brief Copy constructor.
@@ -106,16 +121,26 @@ class Array
        */
       explicit Array( const Array& array );
 
+      /**
+       * \brief Copy constructor with a specific allocator.
+       *
+       * \param array The array to be copied.
+       * \param allocator The allocator to be associated with this array.
+       */
+      explicit Array( const Array& array, const AllocatorType& allocator );
+
       /**
        * \brief Copy constructor.
        *
        * \param array The array to be copied.
        * \param begin The first index which should be copied.
        * \param size The number of elements that should be copied.
+       * \param allocator The allocator to be associated with this array.
        */
       Array( const Array& array,
              IndexType begin,
-             IndexType size = 0 );
+             IndexType size = 0,
+             const AllocatorType& allocator = AllocatorType() );
 
       /**
        * \brief Move constructor for initialization from \e rvalues.
@@ -129,27 +154,40 @@ class Array
        * \ref std::initializer_list, e.g. `{...}`.
        *
        * \param list The initializer list containing elements to be copied.
+       * \param allocator The allocator to be associated with this array.
        */
       template< typename InValue >
-      Array( const std::initializer_list< InValue >& list );
+      Array( const std::initializer_list< InValue >& list,
+             const AllocatorType& allocator = AllocatorType() );
 
       /**
        * \brief Constructor which initializes the array by copying elements from
        * \ref std::list.
        *
        * \param list The STL list containing elements to be copied.
+       * \param allocator The allocator to be associated with this array.
        */
       template< typename InValue >
-      Array( const std::list< InValue >& list );
+      Array( const std::list< InValue >& list,
+             const AllocatorType& allocator = AllocatorType() );
+
 
       /**
        * \brief Constructor which initializes the array by copying elements from
        * \ref std::vector.
        *
        * \param vector The STL vector containing elements to be copied.
+       * \param allocator The allocator to be associated with this array.
        */
       template< typename InValue >
-      Array( const std::vector< InValue >& vector );
+      Array( const std::vector< InValue >& vector,
+             const AllocatorType& allocator = AllocatorType() );
+
+
+      /**
+       * \brief Returns the allocator associated with the array.
+       */
+      AllocatorType getAllocator() const;
 
       /**
        * \brief Returns a \ref String representation of the array type in C++ style.
@@ -586,7 +624,7 @@ class Array
    protected:
 
       /** \brief Method for releasing (deallocating) array data. */
-      void releaseData() const;
+      void releaseData();
 
       /** \brief Number of elements in the array. */
       mutable Index size = 0;
@@ -612,28 +650,33 @@ class Array
        * more arrays. This is to avoid unnecessary dynamic memory allocation.
        */
       mutable int* referenceCounter = nullptr;
+
+      /**
+       * \brief The internal allocator instance.
+       */
+      Allocator allocator;
 };
 
-template< typename Value, typename Device, typename Index >
-std::ostream& operator<<( std::ostream& str, const Array< Value, Device, Index >& array );
+template< typename Value, typename Device, typename Index, typename Allocator >
+std::ostream& operator<<( std::ostream& str, const Array< Value, Device, Index, Allocator >& array );
 
 /**
  * \brief Serialization of arrays into binary files.
  */
-template< typename Value, typename Device, typename Index >
-File& operator<<( File& file, const Array< Value, Device, Index >& array );
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator<<( File& file, const Array< Value, Device, Index, Allocator >& array );
 
-template< typename Value, typename Device, typename Index >
-File& operator<<( File&& file, const Array< Value, Device, Index >& array );
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator<<( File&& file, const Array< Value, Device, Index, Allocator >& array );
 
 /**
  * \brief Deserialization of arrays from binary files.
  */
-template< typename Value, typename Device, typename Index >
-File& operator>>( File& file, Array< Value, Device, Index >& array );
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator>>( File& file, Array< Value, Device, Index, Allocator >& array );
 
-template< typename Value, typename Device, typename Index >
-File& operator>>( File&& file, Array< Value, Device, Index >& array );
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator>>( File&& file, Array< Value, Device, Index, Allocator >& array );
 
 } // namespace Containers
 
diff --git a/src/TNL/Containers/Array.hpp b/src/TNL/Containers/Array.hpp
index 7b4f17f2a38565cd841216ff767e6220d4536664..a667f5a36ca51791f79b65e8eca7e6ef5b7e0997 100644
--- a/src/TNL/Containers/Array.hpp
+++ b/src/TNL/Containers/Array.hpp
@@ -27,39 +27,34 @@ namespace Containers {
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
-Array()
-: size( 0 ),
-  data( 0 ),
-  allocationPointer( 0 ),
-  referenceCounter( 0 )
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
+Array( const Allocator& allocator )
+: allocator( allocator )
 {
 }
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
-Array( const IndexType& size )
-: size( 0 ),
-  data( 0 ),
-  allocationPointer( 0 ),
-  referenceCounter( 0 )
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
+Array( const IndexType& size, const AllocatorType& allocator )
+: allocator( allocator )
 {
    this->setSize( size );
 }
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
 Array( Value* data,
-       const IndexType& size )
-: size( 0 ),
-  data( nullptr ),
-  allocationPointer( nullptr ),
-  referenceCounter( 0 )
+       const IndexType& size,
+       const AllocatorType& allocator )
+: allocator( allocator )
 {
    this->setSize( size );
    Algorithms::ArrayOperations< Device >::copyMemory( this->getData(), data, size );
@@ -67,13 +62,10 @@ Array( Value* data,
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
-Array( const Array< Value, Device, Index >& array )
-: size( 0 ),
-  data( nullptr ),
-  allocationPointer( nullptr ),
-  referenceCounter( 0 )
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
+Array( const Array< Value, Device, Index, Allocator >& array )
 {
    this->setSize( array.getSize() );
    Algorithms::ArrayOperations< Device >::copyMemory( this->getData(), array.getData(), array.getSize() );
@@ -81,15 +73,27 @@ Array( const Array< Value, Device, Index >& array )
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
-Array( const Array< Value, Device, Index >& array,
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
+Array( const Array< Value, Device, Index, Allocator >& array,
+       const AllocatorType& allocator )
+: allocator( allocator )
+{
+   this->setSize( array.getSize() );
+   Algorithms::ArrayOperations< Device >::copyMemory( this->getData(), array.getData(), array.getSize() );
+}
+
+template< typename Value,
+          typename Device,
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
+Array( const Array< Value, Device, Index, Allocator >& array,
        IndexType begin,
-       IndexType size )
-: size( size ),
-  data( nullptr ),
-  allocationPointer( nullptr ),
-  referenceCounter( 0 )
+       IndexType size,
+       const AllocatorType& allocator )
+: allocator( allocator )
 {
    if( size == 0 )
       size = array.getSize() - begin;
@@ -102,14 +106,13 @@ Array( const Array< Value, Device, Index >& array,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename InValue >
-Array< Value, Device, Index >::
-Array( const std::initializer_list< InValue >& list )
-: size( 0 ),
-  data( 0 ),
-  allocationPointer( 0 ),
-  referenceCounter( 0 )
+Array< Value, Device, Index, Allocator >::
+Array( const std::initializer_list< InValue >& list,
+       const AllocatorType& allocator )
+: allocator( allocator )
 {
    this->setSize( list.size() );
    // Here we assume that the underlying array for std::initializer_list is
@@ -120,14 +123,13 @@ Array( const std::initializer_list< InValue >& list )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename InValue >
-Array< Value, Device, Index >::
-Array( const std::list< InValue >& list )
-: size( 0 ),
-  data( 0 ),
-  allocationPointer( 0 ),
-  referenceCounter( 0 )
+Array< Value, Device, Index, Allocator >::
+Array( const std::list< InValue >& list,
+       const AllocatorType& allocator )
+: allocator( allocator )
 {
    this->setSize( list.size() );
    Algorithms::ArrayOperations< Device >::copySTLList( this->getData(), list );
@@ -135,14 +137,13 @@ Array( const std::list< InValue >& list )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename InValue >
-Array< Value, Device, Index >::
-Array( const std::vector< InValue >& vector )
-: size( 0 ),
-  data( 0 ),
-  allocationPointer( 0 ),
-  referenceCounter( 0 )
+Array< Value, Device, Index, Allocator >::
+Array( const std::vector< InValue >& vector,
+       const AllocatorType& allocator )
+: allocator( allocator )
 {
    this->setSize( vector.size() );
    Algorithms::ArrayOperations< Device, Devices::Host >::copyMemory( this->getData(), vector.data(), vector.size() );
@@ -150,9 +151,21 @@ Array( const std::vector< InValue >& vector )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
+Allocator
+Array< Value, Device, Index, Allocator >::
+getAllocator() const
+{
+   return allocator;
+}
+
+template< typename Value,
+          typename Device,
+          typename Index,
+          typename Allocator >
 String
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getType()
 {
    return String( "Containers::Array< " ) +
@@ -163,9 +176,10 @@ getType()
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 String
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getTypeVirtual() const
 {
    return this->getType();
@@ -173,9 +187,10 @@ getTypeVirtual() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 String
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getSerializationType()
 {
    return Algorithms::ArrayIO< Value, Device, Index >::getSerializationType();
@@ -183,9 +198,10 @@ getSerializationType()
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 String
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getSerializationTypeVirtual() const
 {
    return this->getSerializationType();
@@ -193,23 +209,24 @@ getSerializationTypeVirtual() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
-releaseData() const
+Array< Value, Device, Index, Allocator >::
+releaseData()
 {
    if( this->referenceCounter )
    {
       if( --*this->referenceCounter == 0 )
       {
-         Algorithms::ArrayOperations< Device >::freeMemory( this->allocationPointer );
+         allocator.deallocate( this->allocationPointer, this->size );
          delete this->referenceCounter;
          //std::cerr << "Deallocating reference counter " << this->referenceCounter << std::endl;
       }
    }
    else
       if( allocationPointer )
-         Algorithms::ArrayOperations< Device >::freeMemory( this->allocationPointer );
+         allocator.deallocate( this->allocationPointer, this->size );
    this->allocationPointer = 0;
    this->data = 0;
    this->size = 0;
@@ -218,9 +235,10 @@ releaseData() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 setSize( Index size )
 {
    TNL_ASSERT_GE( size, (Index) 0, "Array size must be non-negative." );
@@ -232,7 +250,7 @@ setSize( Index size )
    // Allocating zero bytes is useless. Moreover, the allocators don't behave the same way:
    // "operator new" returns some non-zero address, the latter returns a null pointer.
    if( size > 0 ) {
-      Algorithms::ArrayOperations< Device >::allocateMemory( this->allocationPointer, size );
+      this->allocationPointer = allocator.allocate( size );
       this->data = this->allocationPointer;
       this->size = size;
       TNL_ASSERT_TRUE( this->allocationPointer,
@@ -242,10 +260,11 @@ setSize( Index size )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 Index
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getSize() const
 {
    return this->size;
@@ -253,10 +272,11 @@ getSize() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename ArrayT >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 setLike( const ArrayT& array )
 {
    setSize( array.getSize() );
@@ -264,9 +284,10 @@ setLike( const ArrayT& array )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 bind( Value* data,
       const Index size )
 {
@@ -278,10 +299,11 @@ bind( Value* data,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename ArrayT >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 bind( const ArrayT& array,
       const IndexType& begin,
       const IndexType& size )
@@ -318,10 +340,11 @@ bind( const ArrayT& array,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< int Size >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 bind( StaticArray< Size, Value >& array )
 {
    this->releaseData();
@@ -331,9 +354,10 @@ bind( StaticArray< Size, Value >& array )
 
 template< typename Value,
           typename Device,
-          typename Index >
-typename Array< Value, Device, Index >::ViewType
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+typename Array< Value, Device, Index, Allocator >::ViewType
+Array< Value, Device, Index, Allocator >::
 getView( IndexType begin, IndexType end )
 {
    if( end == 0 )
@@ -343,9 +367,10 @@ getView( IndexType begin, IndexType end )
 
 template< typename Value,
           typename Device,
-          typename Index >
-typename Array< Value, Device, Index >::ConstViewType
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+typename Array< Value, Device, Index, Allocator >::ConstViewType
+Array< Value, Device, Index, Allocator >::
 getView( IndexType begin, IndexType end ) const
 {
    if( end == 0 )
@@ -355,9 +380,10 @@ getView( IndexType begin, IndexType end ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
-typename Array< Value, Device, Index >::ConstViewType
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+typename Array< Value, Device, Index, Allocator >::ConstViewType
+Array< Value, Device, Index, Allocator >::
 getConstView( IndexType begin, IndexType end ) const
 {
    if( end == 0 )
@@ -367,8 +393,9 @@ getConstView( IndexType begin, IndexType end ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
 operator ViewType()
 {
    return getView();
@@ -376,8 +403,9 @@ operator ViewType()
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
 operator ConstViewType() const
 {
    return getConstView();
@@ -385,10 +413,11 @@ operator ConstViewType() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
-swap( Array< Value, Device, Index >& array )
+Array< Value, Device, Index, Allocator >::
+swap( Array< Value, Device, Index, Allocator >& array )
 {
    TNL::swap( this->size, array.size );
    TNL::swap( this->data, array.data );
@@ -398,9 +427,10 @@ swap( Array< Value, Device, Index >& array )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 reset()
 {
    this->releaseData();
@@ -408,10 +438,11 @@ reset()
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 const Value*
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getData() const
 {
    return this->data;
@@ -419,10 +450,11 @@ getData() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 Value*
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getData()
 {
    return this->data;
@@ -430,10 +462,11 @@ getData()
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 const Value*
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getArrayData() const
 {
    return this->data;
@@ -441,10 +474,11 @@ getArrayData() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 Value*
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getArrayData()
 {
    return this->data;
@@ -452,9 +486,10 @@ getArrayData()
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 setElement( const Index& i, const Value& x )
 {
    TNL_ASSERT_GE( i, (Index) 0, "Element index must be non-negative." );
@@ -464,9 +499,10 @@ setElement( const Index& i, const Value& x )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 Value
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 getElement( const Index& i ) const
 {
    TNL_ASSERT_GE( i, (Index) 0, "Element index must be non-negative." );
@@ -476,10 +512,11 @@ getElement( const Index& i ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 inline Value&
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 operator[]( const Index& i )
 {
    TNL_ASSERT_GE( i, (Index) 0, "Element index must be non-negative." );
@@ -489,10 +526,11 @@ operator[]( const Index& i )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 __cuda_callable__
 inline const Value&
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 operator[]( const Index& i ) const
 {
    TNL_ASSERT_GE( i, (Index) 0, "Element index must be non-negative." );
@@ -502,10 +540,11 @@ operator[]( const Index& i ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >&
-Array< Value, Device, Index >::
-operator=( const Array< Value, Device, Index >& array )
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >&
+Array< Value, Device, Index, Allocator >::
+operator=( const Array< Value, Device, Index, Allocator >& array )
 {
    //TNL_ASSERT_EQ( array.getSize(), this->getSize(), "Array sizes must be the same." );
    if( this->getSize() != array.getSize() )
@@ -520,10 +559,11 @@ operator=( const Array< Value, Device, Index >& array )
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >&
-Array< Value, Device, Index >::
-operator=( Array< Value, Device, Index >&& array )
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >&
+Array< Value, Device, Index, Allocator >::
+operator=( Array< Value, Device, Index, Allocator >&& array )
 {
    reset();
 
@@ -542,10 +582,11 @@ operator=( Array< Value, Device, Index >&& array )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename T >
-Array< Value, Device, Index >&
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >&
+Array< Value, Device, Index, Allocator >::
 operator=( const T& data )
 {
    Algorithms::ArrayAssignment< Array, T >::resize( *this, data );
@@ -555,10 +596,11 @@ operator=( const T& data )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename InValue >
-Array< Value, Device, Index >&
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >&
+Array< Value, Device, Index, Allocator >::
 operator=( const std::list< InValue >& list )
 {
    this->setSize( list.size() );
@@ -568,10 +610,11 @@ operator=( const std::list< InValue >& list )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename InValue >
-Array< Value, Device, Index >&
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >&
+Array< Value, Device, Index, Allocator >::
 operator=( const std::vector< InValue >& vector )
 {
    if( this->getSize() != vector.size() )
@@ -582,10 +625,11 @@ operator=( const std::vector< InValue >& vector )
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename ArrayT >
 bool
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 operator==( const ArrayT& array ) const
 {
    if( array.getSize() != this->getSize() )
@@ -600,10 +644,11 @@ operator==( const ArrayT& array ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename ArrayT >
 bool
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 operator!=( const ArrayT& array ) const
 {
    return ! ( *this == array );
@@ -611,9 +656,10 @@ operator!=( const ArrayT& array ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 setValue( const ValueType& v,
           IndexType begin,
           IndexType end )
@@ -626,10 +672,11 @@ setValue( const ValueType& v,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
    template< typename Function >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 evaluate( const Function& f,
           IndexType begin,
           IndexType end )
@@ -640,9 +687,10 @@ evaluate( const Function& f,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 bool
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 containsValue( const ValueType& v,
                IndexType begin,
                IndexType end ) const
@@ -656,9 +704,10 @@ containsValue( const ValueType& v,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 bool
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 containsOnlyValue( const ValueType& v,
                    IndexType begin,
                    IndexType end ) const
@@ -672,10 +721,11 @@ containsOnlyValue( const ValueType& v,
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 bool
 __cuda_callable__
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 empty() const
 {
    return ( data == nullptr );
@@ -683,9 +733,10 @@ empty() const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 save( const String& fileName ) const
 {
    File( fileName, std::ios_base::out ) << *this;
@@ -693,9 +744,10 @@ save( const String& fileName ) const
 
 template< typename Value,
           typename Device,
-          typename Index >
+          typename Index,
+          typename Allocator >
 void
-Array< Value, Device, Index >::
+Array< Value, Device, Index, Allocator >::
 load( const String& fileName )
 {
    File( fileName, std::ios_base::in ) >> *this;
@@ -703,15 +755,16 @@ load( const String& fileName )
 
 template< typename Value,
           typename Device,
-          typename Index >
-Array< Value, Device, Index >::
+          typename Index,
+          typename Allocator >
+Array< Value, Device, Index, Allocator >::
 ~Array()
 {
    this->releaseData();
 }
 
-template< typename Value, typename Device, typename Index >
-std::ostream& operator<<( std::ostream& str, const Array< Value, Device, Index >& array )
+template< typename Value, typename Device, typename Index, typename Allocator >
+std::ostream& operator<<( std::ostream& str, const Array< Value, Device, Index, Allocator >& array )
 {
    str << "[ ";
    if( array.getSize() > 0 )
@@ -725,8 +778,8 @@ std::ostream& operator<<( std::ostream& str, const Array< Value, Device, Index >
 }
 
 // Serialization of arrays into binary files.
-template< typename Value, typename Device, typename Index >
-File& operator<<( File& file, const Array< Value, Device, Index >& array )
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator<<( File& file, const Array< Value, Device, Index, Allocator >& array )
 {
    using IO = Algorithms::ArrayIO< Value, Device, Index >;
    saveObjectType( file, IO::getSerializationType() );
@@ -736,16 +789,16 @@ File& operator<<( File& file, const Array< Value, Device, Index >& array )
    return file;
 }
 
-template< typename Value, typename Device, typename Index >
-File& operator<<( File&& file, const Array< Value, Device, Index >& array )
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator<<( File&& file, const Array< Value, Device, Index, Allocator >& array )
 {
    File& f = file;
    return f << array;
 }
 
 // Deserialization of arrays from binary files.
-template< typename Value, typename Device, typename Index >
-File& operator>>( File& file, Array< Value, Device, Index >& array )
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator>>( File& file, Array< Value, Device, Index, Allocator >& array )
 {
    using IO = Algorithms::ArrayIO< Value, Device, Index >;
    const String type = getObjectType( file );
@@ -760,8 +813,8 @@ File& operator>>( File& file, Array< Value, Device, Index >& array )
    return file;
 }
 
-template< typename Value, typename Device, typename Index >
-File& operator>>( File&& file, Array< Value, Device, Index >& array )
+template< typename Value, typename Device, typename Index, typename Allocator >
+File& operator>>( File&& file, Array< Value, Device, Index, Allocator >& array )
 {
    File& f = file;
    return f >> array;
diff --git a/src/UnitTests/Containers/ArrayTest.h b/src/UnitTests/Containers/ArrayTest.h
index fe0e538b50e9f621f49b051eb6bc010fd7e55eb6..9960e55f337717ab29176f795982f99e2c3d4b7e 100644
--- a/src/UnitTests/Containers/ArrayTest.h
+++ b/src/UnitTests/Containers/ArrayTest.h
@@ -145,7 +145,6 @@ TYPED_TEST_SUITE( ArrayTest, ArrayTypes );
 TYPED_TEST( ArrayTest, constructors )
 {
    using ArrayType = typename TestFixture::ArrayType;
-   //using ArrayType = Array< float, Devices::Host, long >;
 
    ArrayType u;
    EXPECT_EQ( u.getSize(), 0 );
@@ -178,13 +177,61 @@ TYPED_TEST( ArrayTest, constructors )
    EXPECT_EQ( a2.getElement( 2 ), 6 );
 
    std::vector< int > q = { 7, 8, 9 };
-
    ArrayType a3( q );
    EXPECT_EQ( a3.getElement( 0 ), 7 );
    EXPECT_EQ( a3.getElement( 1 ), 8 );
    EXPECT_EQ( a3.getElement( 2 ), 9 );
 }
 
+TYPED_TEST( ArrayTest, constructorsWithAllocators )
+{
+   using ArrayType = typename TestFixture::ArrayType;
+   using AllocatorType = typename ArrayType::AllocatorType;
+
+   AllocatorType allocator;
+
+   ArrayType u( allocator );
+   EXPECT_EQ( u.getAllocator(), allocator );
+   u.setSize( 10 );
+   EXPECT_EQ( u.getSize(), 10 );
+
+   ArrayType v( 10, allocator );
+   EXPECT_EQ( v.getSize(), 10 );
+   EXPECT_EQ( v.getAllocator(), allocator );
+   v.reset();
+   EXPECT_EQ( v.getAllocator(), allocator );
+
+   // deep copy
+   ArrayType w( u, allocator );
+   EXPECT_NE( w.getData(), u.getData() );
+   EXPECT_EQ( w.getSize(), u.getSize() );
+   for( int i = 0; i < 10; i++ )
+      EXPECT_EQ( u.getElement( i ), w.getElement( i ) );
+   EXPECT_EQ( w.getAllocator(), allocator );
+   u.reset();
+   EXPECT_EQ( w.getSize(), 10 );
+
+   ArrayType a1( { 1, 2, 3 }, allocator );
+   EXPECT_EQ( a1.getElement( 0 ), 1 );
+   EXPECT_EQ( a1.getElement( 1 ), 2 );
+   EXPECT_EQ( a1.getElement( 2 ), 3 );
+   EXPECT_EQ( a1.getAllocator(), allocator );
+
+   std::list< int > l = { 4, 5, 6 };
+   ArrayType a2( l, allocator );
+   EXPECT_EQ( a2.getElement( 0 ), 4 );
+   EXPECT_EQ( a2.getElement( 1 ), 5 );
+   EXPECT_EQ( a2.getElement( 2 ), 6 );
+   EXPECT_EQ( a2.getAllocator(), allocator );
+
+   std::vector< int > q = { 7, 8, 9 };
+   ArrayType a3( q, allocator );
+   EXPECT_EQ( a3.getElement( 0 ), 7 );
+   EXPECT_EQ( a3.getElement( 1 ), 8 );
+   EXPECT_EQ( a3.getElement( 2 ), 9 );
+   EXPECT_EQ( a3.getAllocator(), allocator );
+}
+
 TYPED_TEST( ArrayTest, setSize )
 {
    using ArrayType = typename TestFixture::ArrayType;
@@ -519,4 +566,4 @@ TYPED_TEST( ArrayTest, LoadViaView )
 #endif // HAVE_GTEST
 
 
-#include "../main.h"
\ No newline at end of file
+#include "../main.h"