diff --git a/src/TNL/Containers/MultiArray.h b/src/TNL/Containers/MultiArray.h
deleted file mode 100644
index cdcda634bcf768759c44f486790a3511445392d4..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiArray.h
+++ /dev/null
@@ -1,371 +0,0 @@
-/***************************************************************************
-                          MultiArray.h  -  description
-                             -------------------
-    begin                : Nov 25, 2010
-    copyright            : (C) 2010 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-#include <iostream>
-#include <TNL/Containers/Array.h>
-#include <TNL/Containers/StaticVector.h>
-#include <TNL/Assert.h>
-
-namespace TNL {
-namespace Containers {   
-
-template< int Dimension, typename Value = double, typename Device = Devices::Host, typename Index = int >
-class MultiArray : public Array< Value, Device, Index >
-{
-};
-
-template< typename Value, typename Device, typename Index >
-class MultiArray< 1, Value, Device, Index > : public Array< Value, Device, Index >
-{
-   public:
-   enum { Dimension = 1};
-   typedef Value ValueType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiArray< 1, Value, Devices::Host, Index > HostType;
-   typedef MultiArray< 1, Value, Devices::Cuda, Index > CudaType;
-
-
-   MultiArray();
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index iSize );
-
-   void setDimensions( const Containers::StaticVector< 1, Index >& dimensions );
-
-   __cuda_callable__ void getDimensions( Index& iSize ) const;
-
-   __cuda_callable__ const Containers::StaticVector< 1, Index >& getDimensions() const;
-
-   //! Set dimensions of the array using another array as a template
-   template< typename MultiArray >
-   void setLike( const MultiArray& v );
- 
-   void reset();
-
-   __cuda_callable__ Index getElementIndex( const Index i ) const;
-
-   void setElement( const Index i, Value value );
-
-   //! This method can be used for general access to the elements of the arrays.
-   /*! It does not return reference but value. So it can be used to access
-    *  arrays in different address space (usually GPU device).
-    *  See also operator().
-    */
-   Value getElement( const Index i ) const;
-
-   //! Operator for accessing elements of the array.
-   __cuda_callable__ Value& operator()( const Index i );
-
-   __cuda_callable__ const Value& operator()( const Index i ) const;
-
-
-   template< typename MultiArrayT >
-   bool operator == ( const MultiArrayT& array ) const;
-
-   template< typename MultiArrayT >
-   bool operator != ( const MultiArrayT& array ) const;
-
-   MultiArray< 1, Value, Device, Index >& operator = ( const MultiArray< 1, Value, Device, Index >& array );
-
-   template< typename MultiArrayT >
-   MultiArray< 1, Value, Device, Index >& operator = ( const MultiArrayT& array );
-
-   //! Method for saving the object to a file as a binary data
-   bool save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   bool load( File& file );
-
-   bool save( const String& fileName ) const;
-
-   bool load( const String& fileName );
-
-   protected:
-
-   Containers::StaticVector< 1, Index > dimensions;
-};
-
-template< typename Value, typename Device, typename Index >
-class MultiArray< 2, Value, Device, Index > : public Array< Value, Device, Index >
-{
-   public:
-   enum { Dimension = 2 };
-   typedef Value ValueType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiArray< 2, Value, Devices::Host, Index > HostType;
-   typedef MultiArray< 2, Value, Devices::Cuda, Index > CudaType;
-
-
-   MultiArray();
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index jSize, const Index iSize );
-
-   void setDimensions( const Containers::StaticVector< 2, Index >& dimensions );
-
-   __cuda_callable__ void getDimensions( Index& jSize, Index& iSize ) const;
-
-   __cuda_callable__ const Containers::StaticVector< 2, Index >& getDimensions() const;
-
-   //! Set dimensions of the array using another array as a template
-   template< typename MultiArray >
-   void setLike( const MultiArray& v );
-
-   void reset();
-
-   __cuda_callable__ Index getElementIndex( const Index j, const Index i ) const;
-
-   void setElement( const Index j, const Index i, Value value );
-
-   //! This method can be used for general access to the elements of the arrays.
-   /*! It does not return reference but value. So it can be used to access
-    *  arrays in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Value getElement( const Index j, const Index i ) const;
-
-   //! Operator for accessing elements of the array.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of arrays in different address space
-    *  (GPU device usually).
-    */
-   __cuda_callable__ Value& operator()( const Index j, const Index i );
-
-   __cuda_callable__ const Value& operator()( const Index j, const Index i ) const;
-
-   template< typename MultiArrayT >
-   bool operator == ( const MultiArrayT& array ) const;
-
-   template< typename MultiArrayT >
-   bool operator != ( const MultiArrayT& array ) const;
-
-   MultiArray< 2, Value, Device, Index >& operator = ( const MultiArray< 2, Value, Device, Index >& array );
-
-   template< typename MultiArrayT >
-   MultiArray< 2, Value, Device, Index >& operator = ( const MultiArrayT& array );
-
-   //! Method for saving the object to a file as a binary data
-   bool save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   bool load( File& file );
-
-   bool save( const String& fileName ) const;
-
-   bool load( const String& fileName );
-
-   protected:
-
-   Containers::StaticVector< 2, Index > dimensions;
-};
-
-template< typename Value, typename Device, typename Index >
-class MultiArray< 3, Value, Device, Index > : public Array< Value, Device, Index >
-{
-   public:
-
-   enum { Dimension = 3 };
-   typedef Value ValueType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiArray< 3, Value, Devices::Host, Index > HostType;
-   typedef MultiArray< 3, Value, Devices::Cuda, Index > CudaType;
-
-
-   MultiArray();
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index k, const Index j, const Index iSize );
-
-   void setDimensions( const Containers::StaticVector< 3, Index >& dimensions );
-
-   __cuda_callable__ void getDimensions( Index& k, Index& j, Index& iSize ) const;
-
-   __cuda_callable__ const Containers::StaticVector< 3, Index >& getDimensions() const;
-
-   //! Set dimensions of the array using another array as a template
-   template< typename MultiArrayT >
-   void setLike( const MultiArrayT& v );
-
-   void reset();
-
-   __cuda_callable__ Index getElementIndex( const Index k, const Index j, const Index i ) const;
-
-   void setElement( const Index k, const Index j, const Index i, Value value );
-
-   //! This method can be used for general access to the elements of the arrays.
-   /*! It does not return reference but value. So it can be used to access
-    *  arrays in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Value getElement( const Index k, const Index j, const Index i ) const;
-
-   //! Operator for accessing elements of the array.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of arrays in different adress space
-    *  (GPU device usualy).
-    */
-   __cuda_callable__ Value& operator()( const Index k, const Index j, const Index i );
-
-   __cuda_callable__ const Value& operator()( const Index k, const Index j, const Index i ) const;
-
-   template< typename MultiArrayT >
-   bool operator == ( const MultiArrayT& array ) const;
-
-   template< typename MultiArrayT >
-   bool operator != ( const MultiArrayT& array ) const;
-
-   MultiArray< 3, Value, Device, Index >& operator = ( const MultiArray< 3, Value, Device, Index >& array );
-
-   template< typename MultiArrayT >
-   MultiArray< 3, Value, Device, Index >& operator = ( const MultiArrayT& array );
-
-   //! Method for saving the object to a file as a binary data
-   bool save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   bool load( File& file );
-
-   bool save( const String& fileName ) const;
-
-   bool load( const String& fileName );
-
-   protected:
-
-   Containers::StaticVector< 3, Index > dimensions;
-};
-
-template< typename Value, typename Device, typename Index >
-class MultiArray< 4, Value, Device, Index > : public Array< Value, Device, Index >
-{
-   public:
-
-   enum { Dimension = 4 };
-   typedef Value ValueType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiArray< 4, Value, Devices::Host, Index > HostType;
-   typedef MultiArray< 4, Value, Devices::Cuda, Index > CudaType;
-
-
-   MultiArray();
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index l, const Index k, const Index j, const Index iSize );
-
-   void setDimensions( const Containers::StaticVector< 4, Index >& dimensions );
-
-   __cuda_callable__ void getDimensions( Index& l, Index& k, Index& j, Index& iSize ) const;
-
-   __cuda_callable__ const Containers::StaticVector< 4, Index >& getDimensions() const;
-
-   //! Set dimensions of the array using another array as a template
-   template< typename MultiArrayT >
-   void setLike( const MultiArrayT& v );
-
-   void reset();
-
-   __cuda_callable__ Index getElementIndex( const Index l, const Index k, const Index j, const Index i ) const;
-
-   void setElement( const Index l, const Index k, const Index j, const Index i, Value value );
-
-   //! This method can be used for general access to the elements of the arrays.
-   /*! It does not return reference but value. So it can be used to access
-    *  arrays in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Value getElement( const Index l, const Index k, const Index j, const Index i ) const;
-
-   //! Operator for accessing elements of the array.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of arrays in different adress space
-    *  (GPU device usualy).
-    */
-   __cuda_callable__ Value& operator()( const Index l, const Index k, const Index j, const Index i );
-
-   __cuda_callable__ const Value& operator()( const Index l, const Index k, const Index j, const Index i ) const;
-
-   template< typename MultiArrayT >
-   bool operator == ( const MultiArrayT& array ) const;
-
-   template< typename MultiArrayT >
-   bool operator != ( const MultiArrayT& array ) const;
-
-   MultiArray< 4, Value, Device, Index >& operator = ( const MultiArray< 4, Value, Device, Index >& array );
-
-   template< typename MultiArrayT >
-   MultiArray< 4, Value, Device, Index >& operator = ( const MultiArrayT& array );
-
-   //! Method for saving the object to a file as a binary data
-   bool save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   bool load( File& file );
-
-   bool save( const String& fileName ) const;
-
-   bool load( const String& fileName );
-
-   protected:
-
-   Containers::StaticVector< 4, Index > dimensions;
-};
-
-template< typename Value, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 1, Value, device, Index >& array );
-
-template< typename Value, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 2, Value, device, Index >& array );
-
-template< typename Value, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 3, Value, device, Index >& array );
-
-template< typename Value, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 4, Value, device, Index >& array );
-
-} // namespace Containers
-} // namespace TNL
-
-#include <TNL/Containers/MultiArray1D_impl.h>
-#include <TNL/Containers/MultiArray2D_impl.h>
-#include <TNL/Containers/MultiArray3D_impl.h>
-#include <TNL/Containers/MultiArray4D_impl.h>
diff --git a/src/TNL/Containers/MultiArray1D_impl.h b/src/TNL/Containers/MultiArray1D_impl.h
deleted file mode 100644
index 6c8f0b29ceefb37ea4bb237d0a920295603125c0..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiArray1D_impl.h
+++ /dev/null
@@ -1,241 +0,0 @@
-/***************************************************************************
-                          MultiArray1D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once 
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 1, Value, Device, Index > :: MultiArray()
-{
-}
-
-template< typename Value, typename Device, typename Index >
-String MultiArray< 1, Value, Device, Index > :: getType()
-{
-   return String( "Containers::MultiArray< ") +
-          String( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Value >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 1, Value, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 1, Value, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 1, Value, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 1, Value, Device, Index > :: setDimensions( const Index iSize )
-{
-   TNL_ASSERT( iSize > 0,
-              std::cerr << "iSize = " << iSize );
-   dimensions[ 0 ] = iSize;
-   Array< Value, Device, Index >::setSize( iSize );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 1, Value, Device, Index > :: setDimensions( const Containers::StaticVector< 1, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0,
-              std::cerr << " dimensions[ 0 ] = " << dimensions[ 0 ] );
-   this->dimensions = dimensions;
-   Array< Value, Device, Index >::setSize( this->dimensions[ 0 ] );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-void MultiArray< 1, Value, Device, Index > :: setLike( const MultiArrayT& multiArray )
-{
-   setDimensions( multiArray. getDimensions() );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 1, Value, Device, Index >::reset()
-{
-   this->dimensions = Containers::StaticVector< 1, Index >( ( Index ) 0 );
-   Array< Value, Device, Index >::reset();
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-void MultiArray< 1, Value, Device, Index > :: getDimensions( Index& xSize ) const
-{
-   xSize = this->dimensions[ 0 ];
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Containers::StaticVector< 1, Index >& MultiArray< 1, Value, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Index MultiArray< 1, Value, Device, Index > :: getElementIndex( const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ],
-              std::cerr << "i = " << i << " this->dimensions[ 0 ] = " <<  this->dimensions[ 0 ] );
-   return i;
-}
-
-template< typename Value, typename Device, typename Index >
-Value MultiArray< 1, Value, Device, Index > :: getElement( const Index i ) const
-{
-   return Array< Value, Device, Index > :: getElement( getElementIndex( i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 1, Value, Device, Index > :: setElement( const Index i, Value value )
-{
-   Array< Value, Device, Index > :: setElement( getElementIndex( i ), value );
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Value& MultiArray< 1, Value, Device, Index > :: operator()( const Index element )
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( element ) );
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Value& MultiArray< 1, Value, Device, Index > :: operator()( const Index element ) const
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( element ) );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 1, Value, Device, Index > :: operator == ( const MultiArrayT& array ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to compare two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   return Array< Value, Device, Index > :: operator == ( array );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 1, Value, Device, Index > :: operator != ( const MultiArrayT& array ) const
-{
-   return ! ( (* this ) == array );
-}
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 1, Value, Device, Index >&
-   MultiArray< 1, Value, Device, Index > :: operator = ( const MultiArray< 1, Value, Device, Index >& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-MultiArray< 1, Value, Device, Index >&
-   MultiArray< 1, Value, Device, Index > :: operator = ( const MultiArrayT& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 1, Value, Device, Index > :: save( File& file ) const
-{
-   if( ! Array< Value, Device, Index > :: save( file ) )
-   {
-      std::cerr << "I was not able to write the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. save( file ) )
-   {
-      std::cerr << "I was not able to write the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 1, Value, Device, Index > :: load( File& file )
-{
-   if( ! Array< Value, Device, Index > :: load( file ) )
-   {
-      std::cerr << "I was not able to read the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. load( file ) )
-   {
-      std::cerr << "I was not able to read the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 1, Value, Device, Index > :: save( const String& fileName ) const
-{
-   return Object :: save( fileName );
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 1, Value, Device, Index > :: load( const String& fileName )
-{
-   return Object :: load( fileName );
-}
-
-template< typename Value, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 1, Value, Device, Index >& array )
-{
-   for( Index i = 0; i < array. getDimensions()[ 0 ]; i ++ )
-   {
-      str << array. getElement( i ) << " ";
-   }
-   return str;
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiArray2D_impl.h b/src/TNL/Containers/MultiArray2D_impl.h
deleted file mode 100644
index 44d860167968df549a63ed089cad38f9c0919881..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiArray2D_impl.h
+++ /dev/null
@@ -1,255 +0,0 @@
-/***************************************************************************
-                          MultiArray2D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 2, Value, Device, Index > :: MultiArray()
-{
-}
-
-template< typename Value, typename Device, typename Index >
-String MultiArray< 2, Value, Device, Index > :: getType()
-{
-   return String( "Containers::MultiArray< ") +
-          String( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Value >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 2, Value, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 2, Value, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 2, Value, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 2, Value, Device, Index > :: setDimensions( const Index jSize,
-                                                                  const Index iSize )
-{
-   TNL_ASSERT( iSize > 0 && jSize > 0,
-              std::cerr << "iSize = " << iSize
-                   << "jSize = " << jSize );
-
-   dimensions[ 0 ] = iSize;
-   dimensions[ 1 ] = jSize;
-   Array< Value, Device, Index > :: setSize( iSize * jSize );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 2, Value, Device, Index > :: setDimensions( const Containers::StaticVector< 2, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0 && dimensions[ 1 ] > 0,
-              std::cerr << "dimensions = " << dimensions );
-   /****
-    * Swap the dimensions in the tuple to be compatible with the previous method.
-    */
-   this->dimensions. x() = dimensions. y();
-   this->dimensions. y() = dimensions. x();
-   Array< Value, Device, Index > :: setSize( this->dimensions[ 1 ] * this->dimensions[ 0 ] );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-void MultiArray< 2, Value, Device, Index > :: setLike( const MultiArrayT& multiArray )
-{
-   setDimensions( multiArray. getDimensions() );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 2, Value, Device, Index >::reset()
-{
-   this->dimensions = Containers::StaticVector< 2, Index >( ( Index ) 0 );
-   Array< Value, Device, Index >::reset();
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-void MultiArray< 2, Value, Device, Index > :: getDimensions( Index& jSize, Index& iSize ) const
-{
-   iSize = this->dimensions[ 0 ];
-   jSize = this->dimensions[ 1 ];
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Containers::StaticVector< 2, Index >& MultiArray< 2, Value, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Index MultiArray< 2, Value, Device, Index > :: getElementIndex( const Index j, const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ] && j >= 0 && j < this->dimensions[ 1 ],
-              std::cerr << "i = " << i << " j = " << j << " this->dimensions[ 0 ] = " <<  this->dimensions[ 0 ]
-                   << " this->dimensions[ 1 ] = " << this->dimensions[ 1 ] );
-   return j * this->dimensions[ 0 ] + i;
-}
-
-template< typename Value, typename Device, typename Index >
-Value MultiArray< 2, Value, Device, Index > :: getElement( const Index j, const Index i ) const
-{
-   return Array< Value, Device, Index > :: getElement( getElementIndex( j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 2, Value, Device, Index > :: setElement( const Index j, const Index i, Value value )
-{
-   Array< Value, Device, Index > :: setElement( getElementIndex( j, i ), value );
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Value& MultiArray< 2, Value, Device, Index > :: operator()( const Index j, const Index i )
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Value& MultiArray< 2, Value, Device, Index > :: operator()( const Index j, const Index i ) const
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 2, Value, Device, Index > :: operator == ( const MultiArrayT& array ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to compare two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   return Array< Value, Device, Index > :: operator == ( array );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 2, Value, Device, Index > :: operator != ( const MultiArrayT& array ) const
-{
-   return ! ( (* this ) == array );
-}
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 2, Value, Device, Index >&
-   MultiArray< 2, Value, Device, Index > :: operator = ( const MultiArray< 2, Value, Device, Index >& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-MultiArray< 2, Value, Device, Index >&
-   MultiArray< 2, Value, Device, Index > :: operator = ( const MultiArrayT& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 2, Value, Device, Index > :: save( File& file ) const
-{
-   if( ! Array< Value, Device, Index > :: save( file ) )
-   {
-      std::cerr << "I was not able to write the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. save( file ) )
-   {
-      std::cerr << "I was not able to write the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 2, Value, Device, Index > :: load( File& file )
-{
-   if( ! Array< Value, Device, Index > :: load( file ) )
-   {
-      std::cerr << "I was not able to read the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. load( file ) )
-   {
-      std::cerr << "I was not able to read the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 2, Value, Device, Index > :: save( const String& fileName ) const
-{
-   return Object :: save( fileName );
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 2, Value, Device, Index > :: load( const String& fileName )
-{
-   return Object :: load( fileName );
-}
-
-template< typename Value, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 2, Value, Device, Index >& array )
-{
-   for( Index j = 0; j < array. getDimensions()[ 1 ]; j ++ )
-   {
-      for( Index i = 0; i < array. getDimensions()[ 0 ]; i ++ )
-      {
-         str << array. getElement( j, i ) << " ";
-      }
-      str << std::endl;
-   }
-   return str;
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiArray3D_impl.h b/src/TNL/Containers/MultiArray3D_impl.h
deleted file mode 100644
index 9dc3c031795b2e4c5b68fa093512f81911d1abc2..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiArray3D_impl.h
+++ /dev/null
@@ -1,271 +0,0 @@
-/***************************************************************************
-                          MultiArray3D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 3, Value, Device, Index > :: MultiArray()
-{
-}
-
-template< typename Value, typename Device, typename Index >
-String MultiArray< 3, Value, Device, Index > :: getType()
-{
-   return String( "Containers::MultiArray< ") +
-          String( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Value >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 3, Value, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 3, Value, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 3, Value, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 3, Value, Device, Index > :: setDimensions( const Index kSize,
-                                                                       const Index jSize,
-                                                                       const Index iSize )
-{
-   TNL_ASSERT( iSize > 0 && jSize > 0 && kSize > 0,
-              std::cerr << "iSize = " << iSize
-                   << "jSize = " << jSize
-                   << "kSize = " << kSize );
-
-   dimensions[ 0 ] = iSize;
-   dimensions[ 1 ] = jSize;
-   dimensions[ 2 ] = kSize;
-   Array< Value, Device, Index > :: setSize( iSize * jSize * kSize );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 3, Value, Device, Index > :: setDimensions( const Containers::StaticVector< 3, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0 && dimensions[ 1 ] > 0 && dimensions[ 2 ],
-              std::cerr << "dimensions = " << dimensions );
-   /****
-    * Swap the dimensions in the tuple to be compatible with the previous method.
-    */
-   this->dimensions. x() = dimensions. z();
-   this->dimensions. y() = dimensions. y();
-   this->dimensions. z() = dimensions. x();
-   Array< Value, Device, Index > :: setSize( this->dimensions[ 2 ] *
-                                               this->dimensions[ 1 ] *
-                                               this->dimensions[ 0 ] );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-void MultiArray< 3, Value, Device, Index > :: setLike( const MultiArrayT& multiArray )
-{
-   setDimensions( multiArray. getDimensions() );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 3, Value, Device, Index >::reset()
-{
-   this->dimensions = Containers::StaticVector< 3, Index >( ( Index ) 0 );
-   Array< Value, Device, Index >::reset();
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-void MultiArray< 3, Value, Device, Index > :: getDimensions( Index& kSize,
-                                                                  Index& jSize,
-                                                                  Index& iSize ) const
-{
-   iSize = this->dimensions[ 0 ];
-   jSize = this->dimensions[ 1 ];
-   kSize = this->dimensions[ 2 ];
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Containers::StaticVector< 3, Index >& MultiArray< 3, Value, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Index MultiArray< 3, Value, Device, Index > :: getElementIndex( const Index k,
-                                                                     const Index j,
-                                                                     const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ] &&
-              j >= 0 && j < this->dimensions[ 1 ] &&
-              k >= 0 && k < this->dimensions[ 2 ],
-              std::cerr << " i = " << i
-                   << " j = " << j
-                   << " k = " << k
-                   << " this->dimensions = " << this->dimensions );
-   return ( k * this->dimensions[ 1 ]  + j ) * this->dimensions[ 0 ] + i;
-}
-
-template< typename Value, typename Device, typename Index >
-Value MultiArray< 3, Value, Device, Index > :: getElement( const Index k,
-                                                                  const Index j,
-                                                                  const Index i ) const
-{
-   return Array< Value, Device, Index > :: getElement( getElementIndex( k, j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 3, Value, Device, Index > :: setElement( const Index k,
-                                                                    const Index j,
-                                                                    const Index i, Value value )
-{
-   Array< Value, Device, Index > :: setElement( getElementIndex( k, j, i ), value );
-}
-
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Value& MultiArray< 3, Value, Device, Index > :: operator()( const Index k,
-                                                                        const Index j,
-                                                                        const Index i )
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( k, j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Value& MultiArray< 3, Value, Device, Index > :: operator()( const Index k,
-                                                                               const Index j,
-                                                                               const Index i ) const
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( k, j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 3, Value, Device, Index > :: operator == ( const MultiArrayT& array ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to compare two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   return Array< Value, Device, Index > :: operator == ( array );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 3, Value, Device, Index > :: operator != ( const MultiArrayT& array ) const
-{
-   return ! ( (* this ) == array );
-}
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 3, Value, Device, Index >&
-   MultiArray< 3, Value, Device, Index > :: operator = ( const MultiArray< 3, Value, Device, Index >& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-MultiArray< 3, Value, Device, Index >&
-   MultiArray< 3, Value, Device, Index > :: operator = ( const MultiArrayT& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 3, Value, Device, Index > :: save( File& file ) const
-{
-   if( ! Array< Value, Device, Index > :: save( file ) )
-   {
-      std::cerr << "I was not able to write the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. save( file ) )
-   {
-      std::cerr << "I was not able to write the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 3, Value, Device, Index > :: load( File& file )
-{
-   if( ! Array< Value, Device, Index > :: load( file ) )
-   {
-      std::cerr << "I was not able to read the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. load( file ) )
-   {
-      std::cerr << "I was not able to read the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 3, Value, Device, Index >& array )
-{
-   for( Index k = 0; k < array. getDimensions()[ 2 ]; k ++ )
-   {
-      for( Index j = 0; j < array. getDimensions()[ 1 ]; j ++ )
-      {
-         for( Index i = 0; i < array. getDimensions()[ 0 ]; i ++ )
-         {
-            str << array. getElement( k, j, i ) << " ";
-         }
-         str << std::endl;
-      }
-      str << std::endl;
-   }
-   return str;
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiArray4D_impl.h b/src/TNL/Containers/MultiArray4D_impl.h
deleted file mode 100644
index 2b35c1caaf5180645ca1999f69b08a6458cfe4c3..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiArray4D_impl.h
+++ /dev/null
@@ -1,290 +0,0 @@
-/***************************************************************************
-                          MultiArray4D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-   
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 4, Value, Device, Index > :: MultiArray()
-{
-}
-
-template< typename Value, typename Device, typename Index >
-String MultiArray< 4, Value, Device, Index > :: getType()
-{
-   return String( "Containers::MultiArray< ") +
-          String( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Value >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 4, Value, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 4, Value, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Value,
-          typename Device,
-          typename Index >
-String MultiArray< 4, Value, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 4, Value, Device, Index > :: setDimensions( const Index lSize,
-                                                                       const Index kSize,
-                                                                       const Index jSize,
-                                                                       const Index iSize )
-{
-   TNL_ASSERT( iSize > 0 && jSize > 0 && kSize > 0 && lSize > 0,
-              std::cerr << "iSize = " << iSize
-                   << "jSize = " << jSize
-                   << "kSize = " << kSize
-                   << "lSize = " << lSize );
-
-   dimensions[ 0 ] = iSize;
-   dimensions[ 1 ] = jSize;
-   dimensions[ 2 ] = kSize;
-   dimensions[ 3 ] = lSize;
-   Array< Value, Device, Index > :: setSize( iSize * jSize * kSize * lSize );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 4, Value, Device, Index > :: setDimensions( const Containers::StaticVector< 4, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0 && dimensions[ 1 ] > 0 && dimensions[ 2 ] && dimensions[ 3 ] > 0,
-              std::cerr << "dimensions = " << dimensions );
-   /****
-    * Swap the dimensions in the tuple to be compatible with the previous method.
-    */
-   this->dimensions[ 0 ] = dimensions[ 3 ];
-   this->dimensions[ 1 ] = dimensions[ 2 ];
-   this->dimensions[ 2 ] = dimensions[ 1 ];
-   this->dimensions[ 3 ] = dimensions[ 0 ];
-   Array< Value, Device, Index > :: setSize( this->dimensions[ 3 ] *
-                                               this->dimensions[ 2 ] *
-                                               this->dimensions[ 1 ] *
-                                               this->dimensions[ 0 ] );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-void MultiArray< 4, Value, Device, Index > :: setLike( const MultiArrayT& multiArray )
-{
-   setDimensions( multiArray. getDimensions() );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 4, Value, Device, Index >::reset()
-{
-   this->dimensions = Containers::StaticVector< 4, Index >( ( Index ) 0 );
-   Array< Value, Device, Index >::reset();
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-void MultiArray< 4, Value, Device, Index > :: getDimensions( Index& lSize,
-                                                                       Index& kSize,
-                                                                       Index& jSize,
-                                                                       Index& iSize ) const
-{
-   iSize = this->dimensions[ 0 ];
-   jSize = this->dimensions[ 1 ];
-   kSize = this->dimensions[ 2 ];
-   lSize = this->dimensions[ 3 ];
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Containers::StaticVector< 4, Index >& MultiArray< 4, Value, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Index MultiArray< 4, Value, Device, Index > :: getElementIndex( const Index l,
-                                                                     const Index k,
-                                                                     const Index j,
-                                                                     const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ] &&
-              j >= 0 && j < this->dimensions[ 1 ] &&
-              k >= 0 && k < this->dimensions[ 2 ] &&
-              l >= 0 && l < this->dimensions[ 3 ],
-              std::cerr << " i = " << i
-                   << " j = " << j
-                   << " k = " << k
-                   << " l = " << l
-                   << " this->dimensions = " << this->dimensions );
-   return ( ( l * this->dimensions[ 2 ] + k ) * this->dimensions[ 1 ]  + j ) * this->dimensions[ 0 ] + i;
-}
-
-template< typename Value, typename Device, typename Index >
-Value MultiArray< 4, Value, Device, Index > :: getElement( const Index l,
-                                                                       const Index k,
-                                                                       const Index j,
-                                                                       const Index i ) const
-{
-   return Array< Value, Device, Index > :: getElement( getElementIndex( l, k, j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-void MultiArray< 4, Value, Device, Index > :: setElement( const Index l,
-                                                                    const Index k,
-                                                                    const Index j,
-                                                                    const Index i, Value value )
-{
-   Array< Value, Device, Index > :: setElement( getElementIndex( l, k, j, i ), value );
-}
-
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-Value& MultiArray< 4, Value, Device, Index > :: operator()( const Index l,
-                                                                        const Index k,
-                                                                        const Index j,
-                                                                        const Index i )
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( l, k, j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-__cuda_callable__
-const Value& MultiArray< 4, Value, Device, Index > :: operator()( const Index l,
-                                                                               const Index k,
-                                                                               const Index j,
-                                                                               const Index i ) const
-{
-   return Array< Value, Device, Index > :: operator[]( getElementIndex( l, k, j, i ) );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 4, Value, Device, Index > :: operator == ( const MultiArrayT& array ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to compare two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   return Array< Value, Device, Index > :: operator == ( array );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-bool MultiArray< 4, Value, Device, Index > :: operator != ( const MultiArrayT& array ) const
-{
-   return ! ( (* this ) == array );
-}
-
-template< typename Value, typename Device, typename Index >
-MultiArray< 4, Value, Device, Index >&
-   MultiArray< 4, Value, Device, Index > :: operator = ( const MultiArray< 4, Value, Device, Index >& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-   template< typename MultiArrayT >
-MultiArray< 4, Value, Device, Index >&
-   MultiArray< 4, Value, Device, Index > :: operator = ( const MultiArrayT& array )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == array. getDimensions(),
-              std::cerr << "You are attempting to assign two arrays with different dimensions." << std::endl
-                   << "First array dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second array dimensions are ( " << array. getDimensions() << " )" << std::endl; );
-   Array< Value, Device, Index > :: operator = ( array );
-   return ( *this );
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 4, Value, Device, Index > :: save( File& file ) const
-{
-   if( ! Array< Value, Device, Index > :: save( file ) )
-   {
-      std::cerr << "I was not able to write the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. save( file ) )
-   {
-      std::cerr << "I was not able to write the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-bool MultiArray< 4, Value, Device, Index > :: load( File& file )
-{
-   if( ! Array< Value, Device, Index > :: load( file ) )
-   {
-      std::cerr << "I was not able to read the Array of MultiArray." << std::endl;
-      return false;
-   }
-   if( ! dimensions. load( file ) )
-   {
-      std::cerr << "I was not able to read the dimensions of MultiArray." << std::endl;
-      return false;
-   }
-   return true;
-}
-
-template< typename Value, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiArray< 4, Value, Device, Index >& array )
-{
-   for( Index l = 0; l < array. getDimensions()[ 3 ]; l ++ )
-   {
-      for( Index k = 0; k < array. getDimensions()[ 2 ]; k ++ )
-      {
-         for( Index j = 0; j < array. getDimensions()[ 1 ]; j ++ )
-         {
-            for( Index i = 0; i < array. getDimensions()[ 0 ]; i ++ )
-            {
-               str << array. getElement( l, k, j, i ) << " ";
-            }
-            str << std::endl;
-         }
-         str << std::endl;
-      }
-      str << std::endl;
-   }
-   return str;
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiVector.h b/src/TNL/Containers/MultiVector.h
deleted file mode 100644
index eaebf5cd5a3473fb32cb956d3d7ba9fe83edf32a..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiVector.h
+++ /dev/null
@@ -1,369 +0,0 @@
-/***************************************************************************
-                          MultiVector.h  -  description
-                             -------------------
-    begin                : Nov 25, 2010
-    copyright            : (C) 2010 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once 
-
-#include <TNL/Containers/Vector.h>
-#include <TNL/Containers/StaticVector.h>
-#include <TNL/Assert.h>
-
-namespace TNL {
-namespace Containers {   
-   
-template< int Dimension, typename Real = double, typename Device = Devices::Host, typename Index = int >
-class MultiVector : public Vector< Real, Device, Index >
-{
-};
-
-template< typename Real, typename Device, typename Index >
-class MultiVector< 1, Real, Device, Index > : public Vector< Real, Device, Index >
-{
-   public:
-   enum { Dimension = 1};
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiVector< Dimension, Real, Devices::Host, Index > HostType;
-   typedef MultiVector< Dimension, Real, Devices::Cuda, Index > CudaType;
-
-   MultiVector();
-
-   MultiVector( const String& name );
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index iSize );
-
-   void setDimensions( const StaticVector< Dimension, Index >& dimensions );
-
-   void getDimensions( Index& iSize ) const;
-
-   const StaticVector< 1, Index >& getDimensions() const;
-
-   //! Set dimensions of the Vector using another Vector as a template
-   template< typename MultiVector >
-   void setLike( const MultiVector& v );
- 
-   Index getElementIndex( const Index i ) const;
-
-   void setElement( const Index i, Real value );
-
-   //! This method can be used for general access to the elements of the Vectors.
-   /*! It does not return reference but value. So it can be used to access
-    *  Vectors in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Real getElement( const Index i ) const;
-
-   //! Operator for accessing elements of the Vector.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of Vectors in different adress space
-    *  (GPU device usualy).
-    */
-   Real& operator()( const Index i );
-
-   const Real& operator()( const Index i ) const;
-
-   template< typename MultiVector >
-   bool operator == ( const MultiVector& Vector ) const;
-
-   template< typename MultiVector >
-   bool operator != ( const MultiVector& Vector ) const;
-
-   MultiVector< 1, Real, Device, Index >& operator = ( const MultiVector< 1, Real, Device, Index >& Vector );
-
-   template< typename MultiVectorT >
-   MultiVector< 1, Real, Device, Index >& operator = ( const MultiVectorT& Vector );
-
-   //! Method for saving the object to a file as a binary data
-   void save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   void load( File& file );
-
-   void save( const String& fileName ) const;
-
-   void load( const String& fileName );
-
-   protected:
-
-   StaticVector< Dimension, Index > dimensions;
-};
-
-template< typename Real, typename Device, typename Index >
-class MultiVector< 2, Real, Device, Index > : public Vector< Real, Device, Index >
-{
-   public:
-   enum { Dimension = 2 };
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiVector< Dimension, Real, Devices::Host, Index > HostType;
-   typedef MultiVector< Dimension, Real, Devices::Cuda, Index > CudaType;
-
-   MultiVector();
-
-   MultiVector( const String& name );
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index jSize, const Index iSize );
-
-   void setDimensions( const StaticVector< 2, Index >& dimensions );
-
-   void getDimensions( Index& jSize, Index& iSize ) const;
-
-   const StaticVector< 2, Index >& getDimensions() const;
-
-   //! Set dimensions of the Vector using another Vector as a template
-   template< typename MultiVector >
-   void setLike( const MultiVector& v );
-
-   Index getElementIndex( const Index j, const Index i ) const;
-
-   void setElement( const Index j, const Index i, Real value );
-
-   //! This method can be used for general access to the elements of the Vectors.
-   /*! It does not return reference but value. So it can be used to access
-    *  Vectors in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Real getElement( const Index j, const Index i ) const;
-
-   //! Operator for accessing elements of the Vector.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of Vectors in different adress space
-    *  (GPU device usualy).
-    */
-   Real& operator()( const Index j, const Index i );
-
-   const Real& operator()( const Index j, const Index i ) const;
-
-   template< typename MultiVector >
-   bool operator == ( const MultiVector& Vector ) const;
-
-   template< typename MultiVector >
-   bool operator != ( const MultiVector& Vector ) const;
-
-   MultiVector< 2, Real, Device, Index >& operator = ( const MultiVector< 2, Real, Device, Index >& Vector );
-
-   template< typename MultiVectorT >
-   MultiVector< 2, Real, Device, Index >& operator = ( const MultiVectorT& Vector );
-
-   //! Method for saving the object to a file as a binary data
-   void save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   void load( File& file );
-
-   void save( const String& fileName ) const;
-
-   void load( const String& fileName );
-
-   protected:
-
-   StaticVector< 2, Index > dimensions;
-};
-
-template< typename Real, typename Device, typename Index >
-class MultiVector< 3, Real, Device, Index > : public Vector< Real, Device, Index >
-{
-   public:
-
-   enum { Dimension = 3 };
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiVector< Dimension, Real, Devices::Host, Index > HostType;
-   typedef MultiVector< Dimension, Real, Devices::Cuda, Index > CudaType;
-
-   MultiVector();
-
-   MultiVector( const String& name );
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index k, const Index j, const Index iSize );
-
-   void setDimensions( const StaticVector< 3, Index >& dimensions );
-
-   void getDimensions( Index& k, Index& j, Index& iSize ) const;
-
-   const StaticVector< 3, Index >& getDimensions() const;
-
-   //! Set dimensions of the Vector using another Vector as a template
-   template< typename MultiVector >
-   void setLike( const MultiVector& v );
-
-   Index getElementIndex( const Index k, const Index j, const Index i ) const;
-
-   void setElement( const Index k, const Index j, const Index i, Real value );
-
-   //! This method can be used for general access to the elements of the Vectors.
-   /*! It does not return reference but value. So it can be used to access
-    *  Vectors in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Real getElement( const Index k, const Index j, const Index i ) const;
-
-   //! Operator for accessing elements of the Vector.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of Vectors in different adress space
-    *  (GPU device usualy).
-    */
-   Real& operator()( const Index k, const Index j, const Index i );
-
-   const Real& operator()( const Index k, const Index j, const Index i ) const;
-
-   template< typename MultiVector >
-   bool operator == ( const MultiVector& Vector ) const;
-
-   template< typename MultiVector >
-   bool operator != ( const MultiVector& Vector ) const;
-
-   MultiVector< 3, Real, Device, Index >& operator = ( const MultiVector< 3, Real, Device, Index >& Vector );
-
-   template< typename MultiVectorT >
-   MultiVector< 3, Real, Device, Index >& operator = ( const MultiVectorT& Vector );
-
-   //! Method for saving the object to a file as a binary data
-   void save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   void load( File& file );
-
-   void save( const String& fileName ) const;
-
-   void load( const String& fileName );
-
-   protected:
-
-   StaticVector< 3, Index > dimensions;
-};
-
-template< typename Real, typename Device, typename Index >
-class MultiVector< 4, Real, Device, Index > : public Vector< Real, Device, Index >
-{
-   public:
-
-   enum { Dimension = 4 };
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef MultiVector< Dimension, Real, Devices::Host, Index > HostType;
-   typedef MultiVector< Dimension, Real, Devices::Cuda, Index > CudaType;
-
-   MultiVector();
-
-   MultiVector( const String& name );
-
-   static String getType();
-
-   virtual String getTypeVirtual() const;
-
-   static String getSerializationType();
-
-   virtual String getSerializationTypeVirtual() const;
-
-   void setDimensions( const Index l, const Index k, const Index j, const Index iSize );
-
-   void setDimensions( const StaticVector< 4, Index >& dimensions );
-
-   void getDimensions( Index& l, Index& k, Index& j, Index& iSize ) const;
-
-   const StaticVector< 4, Index >& getDimensions() const;
-
-   //! Set dimensions of the Vector using another Vector as a template
-   template< typename MultiVector >
-   void setLike( const MultiVector& v );
-
-   Index getElementIndex( const Index l, const Index k, const Index j, const Index i ) const;
-
-   void setElement( const Index l, const Index k, const Index j, const Index i, Real value );
-
-   //! This method can be used for general access to the elements of the Vectors.
-   /*! It does not return reference but value. So it can be used to access
-    *  Vectors in different adress space (usualy GPU device).
-    *  See also operator().
-    */
-   Real getElement( const Index l, const Index k, const Index j, const Index i ) const;
-
-   //! Operator for accessing elements of the Vector.
-   /*! It returns reference to given elements so it cannot be
-    *  used to access elements of Vectors in different adress space
-    *  (GPU device usualy).
-    */
-   Real& operator()( const Index l, const Index k, const Index j, const Index i );
-
-   const Real& operator()( const Index l, const Index k, const Index j, const Index i ) const;
-
-   template< typename MultiVector >
-   bool operator == ( const MultiVector& Vector ) const;
-
-   template< typename MultiVector >
-   bool operator != ( const MultiVector& Vector ) const;
-
-   MultiVector< 4, Real, Device, Index >& operator = ( const MultiVector< 4, Real, Device, Index >& Vector );
-
-   template< typename MultiVectorT >
-   MultiVector< 4, Real, Device, Index >& operator = ( const MultiVectorT& Vector );
-
-   //! Method for saving the object to a file as a binary data
-   void save( File& file ) const;
-
-   //! Method for restoring the object from a file
-   void load( File& file );
-
-   void save( const String& fileName ) const;
-
-   void load( const String& fileName );
-
-   protected:
-
-   StaticVector< 4, Index > dimensions;
-};
-
-template< typename Real, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 1, Real, device, Index >& Vector );
-
-template< typename Real, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 2, Real, device, Index >& Vector );
-
-template< typename Real, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 3, Real, device, Index >& Vector );
-
-template< typename Real, typename device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 4, Real, device, Index >& Vector );
-
-} // namespace Containers
-} // namespace TNL
-
-#include <TNL/Containers/MultiVector1D_impl.h>
-#include <TNL/Containers/MultiVector2D_impl.h>
-#include <TNL/Containers/MultiVector3D_impl.h>
-#include <TNL/Containers/MultiVector4D_impl.h>
diff --git a/src/TNL/Containers/MultiVector1D_impl.h b/src/TNL/Containers/MultiVector1D_impl.h
deleted file mode 100644
index 1d4f678900a1dafb4213e06b3d0b81d37c44b293..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiVector1D_impl.h
+++ /dev/null
@@ -1,213 +0,0 @@
-/***************************************************************************
-                          MultiVector1D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 1, Real, Device, Index > :: MultiVector()
-{
-}
-
-template< typename Real, typename Device, typename Index >
-String MultiVector< 1, Real, Device, Index > :: getType()
-{
-   return String( "Containers::MultiVector< ") +
-          convertToString( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Real >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 1, Real, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 1, Real, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 1, Real, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: setDimensions( const Index iSize )
-{
-   TNL_ASSERT( iSize > 0,
-              std::cerr << "iSize = " << iSize );
-   dimensions[ 0 ] = iSize;
-   Vector< Real, Device, Index > :: setSize( iSize );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: setDimensions( const StaticVector< Dimension, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0,
-              std::cerr << " dimensions[ 0 ] = " << dimensions[ 0 ] );
-   this->dimensions = dimensions;
-   Vector< Real, Device, Index > :: setSize( this->dimensions[ 0 ] );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-void MultiVector< 1, Real, Device, Index > :: setLike( const MultiVectorT& multiVector )
-{
-   setDimensions( multiVector. getDimensions() );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: getDimensions( Index& xSize ) const
-{
-   xSize = this->dimensions[ 0 ];
-}
-
-template< typename Real, typename Device, typename Index >
-const StaticVector< 1, Index >& MultiVector< 1, Real, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Real, typename Device, typename Index >
-Index MultiVector< 1, Real, Device, Index > :: getElementIndex( const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ],
-              std::cerr << "i = " << i
-                   << "this->dimensions[ 0 ] " << this->dimensions[ 0 ] );
-   return i;
-}
-
-template< typename Real, typename Device, typename Index >
-Real MultiVector< 1, Real, Device, Index > :: getElement( const Index i ) const
-{
-   return Vector< Real, Device, Index > :: getElement( getElementIndex( i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: setElement( const Index i, Real value )
-{
-   Vector< Real, Device, Index > :: setElement( getElementIndex( i ), value );
-}
-
-
-template< typename Real, typename Device, typename Index >
-Real& MultiVector< 1, Real, Device, Index > :: operator()( const Index element )
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( element ) );
-}
-
-template< typename Real, typename Device, typename Index >
-const Real& MultiVector< 1, Real, Device, Index > :: operator()( const Index element ) const
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( element ) );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 1, Real, Device, Index > :: operator == ( const MultiVectorT& vector ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to compare two Vectors with different dimensions." << std::endl
-                   << "First Vector name dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second Vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   return Vector< Real, Device, Index > :: operator == ( vector );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 1, Real, Device, Index > :: operator != ( const MultiVectorT& vector ) const
-{
-   return ! ( (* this ) == vector );
-}
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 1, Real, Device, Index >&
-   MultiVector< 1, Real, Device, Index > :: operator = ( const MultiVector< 1, Real, Device, Index >& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-MultiVector< 1, Real, Device, Index >&
-   MultiVector< 1, Real, Device, Index > :: operator = ( const MultiVectorT& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: save( File& file ) const
-{
-   Vector< Real, Device, Index > :: save( file );
-   dimensions. save( file );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: load( File& file )
-{
-   Vector< Real, Device, Index > :: load( file );
-   dimensions. load( file );
-}
-
-template< typename Real, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 1, Real, Device, Index >& Vector )
-{
-   for( Index i = 0; i < Vector. getDimensions()[ 0 ]; i ++ )
-   {
-      str << Vector. getElement( i ) << " ";
-   }
-   return str;
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: save( const String& fileName ) const
-{
-   Object::save( fileName );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 1, Real, Device, Index > :: load( const String& fileName )
-{
-   Object::load( fileName );
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiVector2D_impl.h b/src/TNL/Containers/MultiVector2D_impl.h
deleted file mode 100644
index a52caff8af18dc470cae3d187a01b8a44a18c877..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiVector2D_impl.h
+++ /dev/null
@@ -1,224 +0,0 @@
-/***************************************************************************
-                          MultiVector2D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 2, Real, Device, Index > :: MultiVector()
-{
-}
-
-template< typename Real, typename Device, typename Index >
-String MultiVector< 2, Real, Device, Index > :: getType()
-{
-   return String( "Containers::MultiVector< ") +
-          convertToString( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Real >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 2, Real, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 2, Real, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 2, Real, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: setDimensions( const Index jSize,
-                                                             const Index iSize )
-{
-   TNL_ASSERT( iSize > 0 && jSize > 0,
-              std::cerr << "iSize = " << iSize
-                   << "jSize = " << jSize );
-
-   dimensions[ 0 ] = iSize;
-   dimensions[ 1 ] = jSize;
-   Vector< Real, Device, Index > :: setSize( iSize * jSize );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: setDimensions( const StaticVector< 2, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0 && dimensions[ 1 ] > 0,
-              std::cerr << "dimensions = " << dimensions );
-   this->dimensions = dimensions;
-   Vector< Real, Device, Index > :: setSize( this->dimensions[ 1 ] * this->dimensions[ 0 ] );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-void MultiVector< 2, Real, Device, Index > :: setLike( const MultiVectorT& multiVector )
-{
-   setDimensions( multiVector. getDimensions() );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: getDimensions( Index& jSize, Index& iSize ) const
-{
-   iSize = this->dimensions[ 0 ];
-   jSize = this->dimensions[ 1 ];
-}
-
-template< typename Real, typename Device, typename Index >
-const StaticVector< 2, Index >& MultiVector< 2, Real, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Real, typename Device, typename Index >
-Index MultiVector< 2, Real, Device, Index > :: getElementIndex( const Index j, const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ] && j >= 0 && j < this->dimensions[ 1 ],
-              std::cerr << "i = " << i
-                   << "j = " << j
-                   << "this->dimensions[ 0 ] = " << this->dimensions[ 0 ]
-                   << "this->dimensions[ 1 ] = " << this->dimensions[ 1 ] );
-   return j * this->dimensions[ 0 ] + i;
-}
-
-template< typename Real, typename Device, typename Index >
-Real MultiVector< 2, Real, Device, Index > :: getElement( const Index j, const Index i ) const
-{
-   return Vector< Real, Device, Index > :: getElement( getElementIndex( j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: setElement( const Index j, const Index i, Real value )
-{
-   Vector< Real, Device, Index > :: setElement( getElementIndex( j, i ), value );
-}
-
-
-template< typename Real, typename Device, typename Index >
-Real& MultiVector< 2, Real, Device, Index > :: operator()( const Index j, const Index i )
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-const Real& MultiVector< 2, Real, Device, Index > :: operator()( const Index j, const Index i ) const
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 2, Real, Device, Index > :: operator == ( const MultiVectorT& vector ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to compare two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   return Vector< Real, Device, Index > :: operator == ( vector );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 2, Real, Device, Index > :: operator != ( const MultiVectorT& vector ) const
-{
-   return ! ( (* this ) == vector );
-}
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 2, Real, Device, Index >&
-   MultiVector< 2, Real, Device, Index > :: operator = ( const MultiVector< 2, Real, Device, Index >& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-MultiVector< 2, Real, Device, Index >&
-   MultiVector< 2, Real, Device, Index > :: operator = ( const MultiVectorT& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: save( File& file ) const
-{
-   Vector< Real, Device, Index > :: save( file );
-   dimensions. save( file );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: load( File& file )
-{
-   Vector< Real, Device, Index > :: load( file );
-   dimensions. load( file );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: save( const String& fileName ) const
-{
-   Object::save( fileName );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 2, Real, Device, Index > :: load( const String& fileName )
-{
-   Object::load( fileName );
-}
-
-template< typename Real, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 2, Real, Device, Index >& Vector )
-{
-   for( Index j = 0; j < Vector. getDimensions()[ 1 ]; j ++ )
-   {
-      for( Index i = 0; i < Vector. getDimensions()[ 0 ]; i ++ )
-      {
-         str << Vector. getElement( j, i ) << " ";
-      }
-      str << std::endl;
-   }
-   return str;
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiVector3D_impl.h b/src/TNL/Containers/MultiVector3D_impl.h
deleted file mode 100644
index 2a280f513434a1ac82ddc54158bf0733877f95c2..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiVector3D_impl.h
+++ /dev/null
@@ -1,248 +0,0 @@
-/***************************************************************************
-                          MultiVector3D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 3, Real, Device, Index > :: MultiVector()
-{
-}
-
-template< typename Real, typename Device, typename Index >
-String MultiVector< 3, Real, Device, Index > :: getType()
-{
-   return String( "Containers::MultiVector< ") +
-          convertToString( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Real >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 3, Real, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 3, Real, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 3, Real, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: setDimensions( const Index kSize,
-                                                             const Index jSize,
-                                                             const Index iSize )
-{
-   TNL_ASSERT( iSize > 0 && jSize > 0 && kSize > 0,
-              std::cerr << "iSize = " << iSize
-                   << "jSize = " << jSize
-                   << "kSize = " << kSize );
-
-   dimensions[ 0 ] = iSize;
-   dimensions[ 1 ] = jSize;
-   dimensions[ 2 ] = kSize;
-   return Vector< Real, Device, Index > :: setSize( iSize * jSize * kSize );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: setDimensions( const StaticVector< 3, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0 && dimensions[ 1 ] > 0 && dimensions[ 2 ],
-              std::cerr << "dimensions = " << dimensions );
-   this->dimensions = dimensions;
-   Vector< Real, Device, Index > :: setSize( this->dimensions[ 2 ] *
-                                             this->dimensions[ 1 ] *
-                                             this->dimensions[ 0 ] );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-void MultiVector< 3, Real, Device, Index > :: setLike( const MultiVectorT& multiVector )
-{
-   setDimensions( multiVector. getDimensions() );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: getDimensions( Index& kSize,
-                                                             Index& jSize,
-                                                             Index& iSize ) const
-{
-   iSize = this->dimensions[ 0 ];
-   jSize = this->dimensions[ 1 ];
-   kSize = this->dimensions[ 2 ];
-}
-
-template< typename Real, typename Device, typename Index >
-const StaticVector< 3, Index >& MultiVector< 3, Real, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Real, typename Device, typename Index >
-Index MultiVector< 3, Real, Device, Index > :: getElementIndex( const Index k,
-                                                                     const Index j,
-                                                                     const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ] &&
-              j >= 0 && j < this->dimensions[ 1 ] &&
-              k >= 0 && k < this->dimensions[ 2 ],
-              std::cerr << " i = " << i
-                   << " j = " << j
-                   << " k = " << k
-                   << " this->dimensions = " << this->dimensions );
-   return ( k * this->dimensions[ 1 ]  + j ) * this->dimensions[ 0 ] + i;
-}
-
-template< typename Real, typename Device, typename Index >
-Real MultiVector< 3, Real, Device, Index > :: getElement( const Index k,
-                                                                  const Index j,
-                                                                  const Index i ) const
-{
-   return Vector< Real, Device, Index > :: getElement( getElementIndex( k, j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: setElement( const Index k,
-                                                                    const Index j,
-                                                                    const Index i, Real value )
-{
-   Vector< Real, Device, Index > :: setElement( getElementIndex( k, j, i ), value );
-}
-
-
-template< typename Real, typename Device, typename Index >
-Real& MultiVector< 3, Real, Device, Index > :: operator()( const Index k,
-                                                                        const Index j,
-                                                                        const Index i )
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( k, j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-const Real& MultiVector< 3, Real, Device, Index > :: operator()( const Index k,
-                                                                               const Index j,
-                                                                               const Index i ) const
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( k, j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 3, Real, Device, Index > :: operator == ( const MultiVectorT& vector ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to compare two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   return Vector< Real, Device, Index > :: operator == ( vector );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 3, Real, Device, Index > :: operator != ( const MultiVectorT& vector ) const
-{
-   return ! ( (* this ) == vector );
-}
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 3, Real, Device, Index >&
-   MultiVector< 3, Real, Device, Index > :: operator = ( const MultiVector< 3, Real, Device, Index >& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-MultiVector< 3, Real, Device, Index >&
-   MultiVector< 3, Real, Device, Index > :: operator = ( const MultiVectorT& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: save( File& file ) const
-{
-   Vector< Real, Device, Index > :: save( file );
-   dimensions. save( file );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: load( File& file )
-{
-   Vector< Real, Device, Index > :: load( file );
-   dimensions. load( file );
-}
-
-template< typename Real, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 3, Real, Device, Index >& Vector )
-{
-   for( Index k = 0; k < Vector. getDimensions()[ 2 ]; k ++ )
-   {
-      for( Index j = 0; j < Vector. getDimensions()[ 1 ]; j ++ )
-      {
-         for( Index i = 0; i < Vector. getDimensions()[ 0 ]; i ++ )
-         {
-            str << Vector. getElement( k, j, i ) << " ";
-         }
-         str << std::endl;
-      }
-      str << std::endl;
-   }
-   return str;
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: save( const String& fileName ) const
-{
-   Object::save( fileName );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 3, Real, Device, Index > :: load( const String& fileName )
-{
-   Object::load( fileName );
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/TNL/Containers/MultiVector4D_impl.h b/src/TNL/Containers/MultiVector4D_impl.h
deleted file mode 100644
index f3c838155da1525fca2304b167c7702abaffe496..0000000000000000000000000000000000000000
--- a/src/TNL/Containers/MultiVector4D_impl.h
+++ /dev/null
@@ -1,269 +0,0 @@
-/***************************************************************************
-                          MultiVector4D_impl.h  -  description
-                             -------------------
-    begin                : Nov 13, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-namespace TNL {
-namespace Containers {   
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 4, Real, Device, Index > :: MultiVector()
-{
-}
-
-template< typename Real, typename Device, typename Index >
-String MultiVector< 4, Real, Device, Index > :: getType()
-{
-   return String( "Containers::MultiVector< ") +
-          convertToString( Dimension ) +
-          String( ", " ) +
-          String( TNL::getType< Real >() ) +
-          String( ", " ) +
-          String( Device :: getDeviceType() ) +
-          String( ", " ) +
-          String( TNL::getType< Index >() ) +
-          String( " >" );
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 4, Real, Device, Index > :: getTypeVirtual() const
-{
-   return this->getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 4, Real, Device, Index > :: getSerializationType()
-{
-   return HostType::getType();
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-String MultiVector< 4, Real, Device, Index > :: getSerializationTypeVirtual() const
-{
-   return this->getSerializationType();
-};
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: setDimensions( const Index lSize,
-                                                             const Index kSize,
-                                                             const Index jSize,
-                                                             const Index iSize )
-{
-   TNL_ASSERT( iSize > 0 && jSize > 0 && kSize > 0 && lSize > 0,
-              std::cerr << "iSize = " << iSize
-                   << "jSize = " << jSize
-                   << "kSize = " << kSize
-                   << "lSize = " << lSize );
-
-   dimensions[ 0 ] = iSize;
-   dimensions[ 1 ] = jSize;
-   dimensions[ 2 ] = kSize;
-   dimensions[ 3 ] = lSize;
-   Vector< Real, Device, Index > :: setSize( iSize * jSize * kSize * lSize );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: setDimensions( const StaticVector< 4, Index >& dimensions )
-{
-   TNL_ASSERT( dimensions[ 0 ] > 0 && dimensions[ 1 ] > 0 && dimensions[ 2 ] && dimensions[ 3 ] > 0,
-              std::cerr << "dimensions = " << dimensions );
-   this->dimensions = dimensions;
-   Vector< Real, Device, Index > :: setSize( this->dimensions[ 3 ] *
-                                             this->dimensions[ 2 ] *
-                                             this->dimensions[ 1 ] *
-                                             this->dimensions[ 0 ] );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-void MultiVector< 4, Real, Device, Index > :: setLike( const MultiVectorT& multiVector )
-{
-   setDimensions( multiVector. getDimensions() );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: getDimensions( Index& lSize,
-                                                                       Index& kSize,
-                                                                       Index& jSize,
-                                                                       Index& iSize ) const
-{
-   iSize = this->dimensions[ 0 ];
-   jSize = this->dimensions[ 1 ];
-   kSize = this->dimensions[ 2 ];
-   lSize = this->dimensions[ 3 ];
-}
-
-template< typename Real, typename Device, typename Index >
-const StaticVector< 4, Index >& MultiVector< 4, Real, Device, Index > :: getDimensions() const
-{
-   return this->dimensions;
-}
-
-template< typename Real, typename Device, typename Index >
-Index MultiVector< 4, Real, Device, Index > :: getElementIndex( const Index l,
-                                                                          const Index k,
-                                                                          const Index j,
-                                                                          const Index i ) const
-{
-   TNL_ASSERT( i >= 0 && i < this->dimensions[ 0 ] &&
-              j >= 0 && j < this->dimensions[ 1 ] &&
-              k >= 0 && k < this->dimensions[ 2 ] &&
-              l >= 0 && l < this->dimensions[ 3 ],
-              std::cerr << " i = " << i
-                   << " j = " << j
-                   << " k = " << k
-                   << " l = " << l
-                   << " this->dimensions = " << this->dimensions );
-   return ( ( l * this->dimensions[ 2 ] + k ) * this->dimensions[ 1 ]  + j ) * this->dimensions[ 0 ] + i;
-}
-
-template< typename Real, typename Device, typename Index >
-Real
-MultiVector< 4, Real, Device, Index >::
-getElement( const Index l,
-            const Index k,
-            const Index j,
-            const Index i ) const
-{
-   return Vector< Real, Device, Index > :: getElement( getElementIndex( l, k, j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: setElement( const Index l,
-                                                                    const Index k,
-                                                                    const Index j,
-                                                                    const Index i, Real value )
-{
-   Vector< Real, Device, Index > :: setElement( getElementIndex( l, k, j, i ), value );
-}
-
-
-template< typename Real, typename Device, typename Index >
-Real&
-MultiVector< 4, Real, Device, Index >::
-operator()( const Index l,
-            const Index k,
-            const Index j,
-            const Index i )
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( l, k, j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-const Real& MultiVector< 4, Real, Device, Index > :: operator()( const Index l,
-                                                                               const Index k,
-                                                                               const Index j,
-                                                                               const Index i ) const
-{
-   return Vector< Real, Device, Index > :: operator[]( getElementIndex( l, k, j, i ) );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 4, Real, Device, Index > :: operator == ( const MultiVectorT& vector ) const
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to compare two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   return Vector< Real, Device, Index > :: operator == ( vector );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-bool MultiVector< 4, Real, Device, Index > :: operator != ( const MultiVectorT& vector ) const
-{
-   return ! ( (* this ) == vector );
-}
-
-template< typename Real, typename Device, typename Index >
-MultiVector< 4, Real, Device, Index >&
-   MultiVector< 4, Real, Device, Index > :: operator = ( const MultiVector< 4, Real, Device, Index >& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-   template< typename MultiVectorT >
-MultiVector< 4, Real, Device, Index >&
-   MultiVector< 4, Real, Device, Index > :: operator = ( const MultiVectorT& vector )
-{
-   // TODO: Static assert on dimensions
-   TNL_ASSERT( this->getDimensions() == vector. getDimensions(),
-              std::cerr << "You are attempting to assign two Vectors with different dimensions." << std::endl
-                   << "First vector dimensions are ( " << this->getDimensions() << " )" << std::endl
-                   << "Second vector dimensions are ( " << vector. getDimensions() << " )" << std::endl; );
-   Vector< Real, Device, Index > :: operator = ( vector );
-   return ( *this );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: save( File& file ) const
-{
-   Vector< Real, Device, Index > :: save( file );
-   dimensions. save( file );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: load( File& file )
-{
-   Vector< Real, Device, Index > :: load( file );
-   dimensions. load( file );
-}
-
-template< typename Real, typename Device, typename Index >
-std::ostream& operator << ( std::ostream& str, const MultiVector< 4, Real, Device, Index >& Vector )
-{
-   for( Index l = 0; l < Vector. getDimensions()[ 3 ]; l ++ )
-   {
-      for( Index k = 0; k < Vector. getDimensions()[ 2 ]; k ++ )
-      {
-         for( Index j = 0; j < Vector. getDimensions()[ 1 ]; j ++ )
-         {
-            for( Index i = 0; i < Vector. getDimensions()[ 0 ]; i ++ )
-            {
-               str << Vector. getElement( l, k, j, i ) << " ";
-            }
-            str << std::endl;
-         }
-         str << std::endl;
-      }
-      str << std::endl;
-   }
-   return str;
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: save( const String& fileName ) const
-{
-   Object::save( fileName );
-}
-
-template< typename Real, typename Device, typename Index >
-void MultiVector< 4, Real, Device, Index > :: load( const String& fileName )
-{
-   Object::load( fileName );
-}
-
-} // namespace Containers
-} // namespace TNL
diff --git a/src/Tools/tnl-diff.h b/src/Tools/tnl-diff.h
index 364dbd294e6fafdd0f417af6142cad2b6f01558d..01be8959fd31f4e8801042738cce995527bc4715 100644
--- a/src/Tools/tnl-diff.h
+++ b/src/Tools/tnl-diff.h
@@ -452,10 +452,6 @@ bool setIndexType( const MeshPointer& meshPointer,
                    const Config::ParameterContainer& parameters )
 {
    String indexType;
-   if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-       parsedObjectType[ 0 ] == "tnlMultiVector" ||                       // TODO: remove deprecated type names
-       parsedObjectType[ 0 ] == "tnlSharedMultiVector"   )                //
-      indexType = parsedObjectType[ 4 ];
    if( parsedObjectType[ 0 ] == "Containers::Vector" ||
        parsedObjectType[ 0 ] == "tnlSharedVector" ||                     // TODO: remove deprecated type names
        parsedObjectType[ 0 ] == "tnlVector" )                            //
@@ -532,10 +528,6 @@ bool setValueType( const MeshPointer& meshPointer,
 {
    String elementType;
 
-   if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-       parsedObjectType[ 0 ] == "tnlMultiVector" ||                         // TODO: remove deprecated type names
-       parsedObjectType[ 0 ] == "tnlSharedMultiVector" )                    //
-      elementType = parsedObjectType[ 2 ];
    if( parsedObjectType[ 0 ] == "Functions::MeshFunction" ||
        parsedObjectType[ 0 ] == "tnlMeshFunction" )                         // TODO: remove deprecated type names
       elementType = parsedObjectType[ 3 ];
diff --git a/src/Tools/tnl-view.h b/src/Tools/tnl-view.h
index cfd958aea319486ea07dcf85cddd661f4e1256cc..d7a25fd1ed67b422a13ab9a13c79a987521ee3e4 100644
--- a/src/Tools/tnl-view.h
+++ b/src/Tools/tnl-view.h
@@ -16,7 +16,6 @@
 #include <TNL/Config/ParameterContainer.h>
 #include <TNL/String.h>
 #include <TNL/Containers/Vector.h>
-#include <TNL/Containers/MultiVector.h>
 #include <TNL/Meshes/Grid.h>
 #include <TNL/Functions/MeshFunction.h>
 #include <TNL/Functions/VectorField.h>
@@ -264,22 +263,6 @@ bool convertObject( const MeshPointer& meshPointer,
       mf.bind( meshPointer, vector );
       mf.write( outputFileName, outputFormat );
    }
-
-   if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-       parsedObjectType[ 0 ] == "tnlMultiVector" ||      // TODO: remove deprecated type names  
-       parsedObjectType[ 0 ] == "tnlSharedMultiVector" ) //
-   {
-      Containers::MultiVector< Dimension, Value, Devices::Host, Index > multiVector;
-      multiVector. load( inputFileName );
-      typedef Meshes::Grid< Dimension, Real, Devices::Host, Index > GridType;
-      typedef typename GridType::PointType PointType;
-      typedef typename GridType::CoordinatesType CoordinatesType;
-//      GridType grid;
-//      grid. setDomain( PointType( 0.0 ), PointType( 1.0 ) );
-//      grid. setDimensions( CoordinatesType( multiVector. getDimensions() ) );
-//      if( ! grid. write( multiVector, outputFileName, outputFormat ) )
-//         return false;
-   }
    return true;
 }
 
@@ -290,10 +273,6 @@ bool setDimension( const MeshPointer& meshPointer,
                     const Config::ParameterContainer& parameters )
 {
    int dimensions( 0 );
-   if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-       parsedObjectType[ 0 ] == "tnlMultiVector" ||                     // TODO: remove deprecated type names
-       parsedObjectType[ 0 ] == "tnlSharedMultiVector" )                //
-      dimensions = atoi( parsedObjectType[ 1 ]. getString() );
    if( parsedObjectType[ 0 ] == "Containers::Vector" ||
        parsedObjectType[ 0 ] == "tnlVector" ||                          // TODO: remove deprecated type names
        parsedObjectType[ 0 ] == "tnlSharedVector" )                     //
@@ -318,10 +297,6 @@ bool setIndexType( const MeshPointer& meshPointer,
                    const Config::ParameterContainer& parameters )
 {
    String indexType;
-   if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-       parsedObjectType[ 0 ] == "tnlMultiVector" ||                        // TODO: remove deprecated type names
-       parsedObjectType[ 0 ] == "tnlSharedMultiVector" )                   //
-      indexType = parsedObjectType[ 4 ];
    if( parsedObjectType[ 0 ] == "Containers::Vector" || 
        parsedObjectType[ 0 ] == "tnlSharedVector" ||                       // TODO: remove deprecated type names
        parsedObjectType[ 0 ] == "tnlVector" )                              //
@@ -395,10 +370,6 @@ bool setValueType( const MeshPointer& meshPointer,
    String elementType;
 
    // TODO: Fix this even for arrays
-   if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-       parsedObjectType[ 0 ] == "tnlMultiVector" ||                            // TODO: remove deprecated type names
-       parsedObjectType[ 0 ] == "tnlSharedMultiVector" )                       //
-      elementType = parsedObjectType[ 2 ];
    if( parsedObjectType[ 0 ] == "Containers::Vector" ||
        parsedObjectType[ 0 ] == "tnlSharedVector" ||                           // TODO: remove deprecated type names
        parsedObjectType[ 0 ] == "tnlVector" )                                  //
@@ -497,11 +468,8 @@ struct FilesProcessor
             error = true;
             continue;
          }
-         if( parsedObjectType[ 0 ] == "Containers::MultiVector" ||
-             parsedObjectType[ 0 ] == "Containers::Vector" ||
-             parsedObjectType[ 0 ] == "tnlMultiVector" ||                     // TODO: remove deprecated type names
-             parsedObjectType[ 0 ] == "tnlSharedMultiVector" ||               // 
-             parsedObjectType[ 0 ] == "tnlSharedVector" ||                    //
+         if( parsedObjectType[ 0 ] == "Containers::Vector" ||
+             parsedObjectType[ 0 ] == "tnlSharedVector" ||                    // TODO: remove deprecated type names
              parsedObjectType[ 0 ] == "tnlVector" )                           //
             setValueType< MeshPointer >( meshPointer, inputFiles[ i ], parsedObjectType, parameters );
          if( parsedObjectType[ 0 ] == "Functions::MeshFunction" ||
diff --git a/src/UnitTests/Containers/CMakeLists.txt b/src/UnitTests/Containers/CMakeLists.txt
index fe75ed458bad9e8e8314a4e351e86dcfc6870bea..99e4d677bc3df482484bb56b87cea9c71c3a6f15 100644
--- a/src/UnitTests/Containers/CMakeLists.txt
+++ b/src/UnitTests/Containers/CMakeLists.txt
@@ -58,16 +58,6 @@ ADD_EXECUTABLE( StaticVectorTest StaticVectorTest.cpp )
 TARGET_COMPILE_OPTIONS( StaticVectorTest PRIVATE ${CXX_TESTS_FLAGS} )
 TARGET_LINK_LIBRARIES( StaticVectorTest ${GTEST_BOTH_LIBRARIES} )
 
-#IF( BUILD_CUDA )
-#   CUDA_ADD_EXECUTABLE( MultiArrayTest MultiArrayTest.cu
-#                        OPTIONS ${CXX_TESTS_FLAGS} )
-#   TARGET_LINK_LIBRARIES( MultiArrayTest ${GTEST_BOTH_LIBRARIES} )
-#ELSE(  BUILD_CUDA )
-#   ADD_EXECUTABLE( MultiArrayTest MultiArrayTest.cpp )
-#   TARGET_COMPILE_OPTIONS( MultiArrayTest PRIVATE ${CXX_TESTS_FLAGS} )
-#   TARGET_LINK_LIBRARIES( MultiArrayTest ${GTEST_BOTH_LIBRARIES} )
-#ENDIF( BUILD_CUDA )
-
 
 ADD_TEST( ListTest ${EXECUTABLE_OUTPUT_PATH}/ListTest${CMAKE_EXECUTABLE_SUFFIX} )
 ADD_TEST( ArrayOperationsTest ${EXECUTABLE_OUTPUT_PATH}/ArrayOperationsTest${CMAKE_EXECUTABLE_SUFFIX} )
@@ -82,7 +72,6 @@ ENDIF()
 ADD_TEST( MultireductionTest ${EXECUTABLE_OUTPUT_PATH}/MultireductionTest${CMAKE_EXECUTABLE_SUFFIX} )
 ADD_TEST( StaticArrayTest ${EXECUTABLE_OUTPUT_PATH}/StaticArrayTest${CMAKE_EXECUTABLE_SUFFIX} )
 ADD_TEST( StaticVectorTest ${EXECUTABLE_OUTPUT_PATH}/StaticVectorTest${CMAKE_EXECUTABLE_SUFFIX} )
-#ADD_TEST( MultiArrayTest ${EXECUTABLE_OUTPUT_PATH}/MultiArrayTest${CMAKE_EXECUTABLE_SUFFIX} )
 
 
 ADD_SUBDIRECTORY( Multimaps )
diff --git a/src/UnitTests/Containers/MultiArrayTest.cpp b/src/UnitTests/Containers/MultiArrayTest.cpp
deleted file mode 100644
index 18e7453ccf928f288127d2de469f60d63005ca7b..0000000000000000000000000000000000000000
--- a/src/UnitTests/Containers/MultiArrayTest.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-/***************************************************************************
-                          MultiArrayTest.cpp  -  description
-                             -------------------
-    begin                : Feb 3, 2014
-    copyright            : (C) 2014 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#include "MultiArrayTest.h"
\ No newline at end of file
diff --git a/src/UnitTests/Containers/MultiArrayTest.cu b/src/UnitTests/Containers/MultiArrayTest.cu
deleted file mode 100644
index 97d7ff312b645d87bbaa230cb355a883678508d1..0000000000000000000000000000000000000000
--- a/src/UnitTests/Containers/MultiArrayTest.cu
+++ /dev/null
@@ -1,11 +0,0 @@
-/***************************************************************************
-                          MultiArrayTest.cu  -  description
-                             -------------------
-    begin                : Feb 3, 2014
-    copyright            : (C) 2014 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#include "MultiArrayTest.h"
diff --git a/src/UnitTests/Containers/MultiArrayTest.h b/src/UnitTests/Containers/MultiArrayTest.h
deleted file mode 100644
index 3e8d330a1cc3f0e2eb7e0539485ba1b737591595..0000000000000000000000000000000000000000
--- a/src/UnitTests/Containers/MultiArrayTest.h
+++ /dev/null
@@ -1,237 +0,0 @@
-/***************************************************************************
-                          MultiArrayTester.h -  description
-                             -------------------
-    begin                : Jul 4, 2012
-    copyright            : (C) 2012 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-#include <TNL/Containers/MultiArray.h>
-
-#ifdef HAVE_GTEST 
-#include "gtest/gtest.h"
-#endif
-
-using namespace TNL;
-using namespace TNL::Containers;
-
-#ifdef HAVE_CUDA
-template< typename ValueType, typename IndexType >
-__global__ void testSetGetElementKernel( MultiArray< 1, ValueType, Devices::Cuda, IndexType >* u )
-{
-   if( threadIdx.x < ( *u ).getDimensions().x() )
-      ( *u )( threadIdx.x ) = threadIdx.x;
-}
-
-template< typename ValueType, typename IndexType >
-__global__ void testSetGetElementKernel( MultiArray< 2, ValueType, Devices::Cuda, IndexType >* u )
-{
-   if( threadIdx.x < ( *u ).getDimensions().x() &&
-       threadIdx.x < ( *u ).getDimensions().y() )
-      ( *u )( threadIdx.x, threadIdx.x ) = threadIdx.x;
-}
-
-template< typename ValueType, typename IndexType >
-__global__ void testSetGetElementKernel( MultiArray< 3, ValueType, Devices::Cuda, IndexType >* u )
-{
-   if( threadIdx.x < ( *u ).getDimensions().x() &&
-       threadIdx.x < ( *u ).getDimensions().y() &&
-       threadIdx.x < ( *u ).getDimensions().z() )
-      ( *u )( threadIdx.x, threadIdx.x, threadIdx.x ) = threadIdx.x;
-}
-
-#endif /* HAVE_CUDA */
-
-#ifdef HAVE_GTEST
-
-TEST( MultiArrayTest, testConstructorDestructor )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u;
-}
-
-TEST( MultiArrayTest, testSetSize )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u, v;
-   u. setDimensions( 10 );
-   v. setDimensions( 10 );
-}
-
-void setDiagonalElement( Containers::MultiArray< 1, ValueType, Device, IndexType >& u,
-                         const IndexType& i,
-                         const ValueType& v )
-{
-   u.setElement( i, v );
-}
-
-void setDiagonalElement( Containers::MultiArray< 2, ValueType, Device, IndexType >& u,
-                         const IndexType& i,
-                         const ValueType& v )
-{
-   u.setElement( i, i, v );
-}
-
-void setDiagonalElement( Containers::MultiArray< 3, ValueType, Device, IndexType >& u,
-                         const IndexType& i,
-                         const ValueType& v )
-{
-   u.setElement( i, i, i, v );
-}
-
-IndexType getDiagonalElement( Containers::MultiArray< 1, ValueType, Device, IndexType >& u,
-                              const IndexType& i )
-{
-   return u.getElement( i );
-}
-
-IndexType getDiagonalElement( Containers::MultiArray< 2, ValueType, Device, IndexType >& u,
-                              const IndexType& i )
-{
-   return u.getElement( i, i );
-}
-
-IndexType getDiagonalElement( Containers::MultiArray< 3, ValueType, Device, IndexType >& u,
-                              const IndexType& i )
-{
-   return u.getElement( i, i, i );
-}
-
-
-TEST( MultiArrayTest, testSetGetElement )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u;
-   u. setDimensions( 10 );
-   if( std::is_same< Device, Devices::Host >::value )
-   {
-      for( int i = 0; i < 10; i ++ )
-         this->setDiagonalElement( u, i, i  );
-   }
-   if( std::is_same< Device, Devices::Cuda >::value )
-   {
-#ifdef HAVE_CUDA
-      MultiArray< Dimension, ValueType, Device, IndexType >* kernel_u =
-               Devices::Cuda::passToDevice( u );
-      testSetGetElementKernel<<< 1, 16 >>>( kernel_u );
-      Devices::Cuda::freeFromDevice( kernel_u );
-      ASSERT_TRUE( TNL_CHECK_CUDA_DEVICE );
-#endif
-   }
-   for( int i = 0; i < 10; i ++ )
-      ASSERT_EQ( getDiagonalElement( u, i ), i );
-};
-
-TEST( MultiArrayTest, testComparisonOperator )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u, v, w;
-   u.setDimensions( 10 );
-   v.setDimensions( 10 );
-   w.setDimensions( 10 );
-   u.setValue( 0 );
-   v.setValue( 0 );
-   w.setValue( 0 );
-   for( int i = 0; i < 10; i ++ )
-   {
-      setDiagonalElement( u, i, i );
-      setDiagonalElement( v, i, i );
-      setDiagonalElement( w, i, 2*1 );
-   }
-   ASSERT_TRUE( u == v );
-   ASSERT_FALSE( u != v );
-   ASSERT_TRUE( u != w );
-   ASSERT_FALSE( u == w );
-};
-
-TEST( MultiArrayTest, testEquivalenceOperator )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u;
-   MultiArray< Dimension, ValueType, Device, IndexType > v;
-   u. setDimensions( 10 );
-   v. setDimensions( 10 );
-   for( int i = 0; i < 10; i ++ )
-      setDiagonalElement( u, i, i );
-   v = u;
-   ASSERT_TRUE( u == v );
-   ASSERT_FALSE( u != v );
-};
-
-TEST( MultiArrayTest, testGetSize )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u;
-   const int maxSize = 10;
-   for( int i = 1; i < maxSize; i ++ )
-      u. setDimensions( i );
-
-   ASSERT_EQ( u. getDimensions().x(), maxSize - 1 );
-};
-
-TEST( MultiArrayTest, testReset )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > u;
-   u.setDimensions( 100 );
-   ASSERT_EQ( u. getDimensions().x(), 100 );
-   u.reset();
-   ASSERT_EQ( u. getDimensions().x(), 0 );
-   u.setDimensions( 100 );
-   ASSERT_EQ( u. getDimensions().x(), 100 );
-   u.reset();
-   ASSERT_EQ( u. getDimensions().x(), 0 );
-
-};
-
-TEST( MultiArrayTest, testSetSizeAndDestructor )
-{
-   using namespace TNL::Containers;
-   for( int i = 1; i < 100; i ++ )
-   {
-      MultiArray< Dimension, ValueType, Device, IndexType > u;
-      u. setDimensions( i );
-   }
-}
-
-TEST( MultiArrayTest, testSaveAndLoad )
-{
-   using namespace TNL::Containers;
-   MultiArray< Dimension, ValueType, Device, IndexType > v;
-   const int size( 10 );
-   ASSERT_TRUE( v. setDimensions( size ) );
-   for( int i = 0; i < size; i ++ )
-      setDiagonalElement( v, i, 3.14147 );
-   File file;
-   file. open( "test-file.tnl", File::Mode::Out );
-   ASSERT_TRUE( v. save( file ) );
-   file. close();
-   MultiArray< Dimension, ValueType, Device, IndexType > u;
-   file. open( "test-file.tnl", File::Mode::In );
-   ASSERT_TRUE( u. load( file ) );
-   file. close();
-   ASSERT_TRUE( u == v );
-
-   EXPECT_EQ( std::remove( "test-file.tnl" ), 0 );
-}
-#endif /* HAVE_GTEST */
-
-int main( int argc, char* argv[] )
-{
-#ifdef HAVE_GTEST
-   ::testing::InitGoogleTest( &argc, argv );
-   return RUN_ALL_TESTS();
-#else
-   return EXIT_FAILURE;
-#endif
-}
-
-
-
-
-
-