table.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_table_h
#define dealii_table_h

#include <deal.II/base/config.h>

#include <deal.II/base/aligned_vector.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/linear_index_iterator.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table_indices.h>

#include <algorithm>
#include <cstddef>

DEAL_II_NAMESPACE_OPEN

// forward declaration
template <int N, typename T>
class TableBase;
template <int N, typename T>
class Table;
template <typename T>
class TransposeTable;
template <typename T>
class Table<1, T>;
template <typename T>
class Table<2, T>;
template <typename T>
class Table<3, T>;
template <typename T>
class Table<4, T>;
template <typename T>
class Table<5, T>;
template <typename T>
class Table<6, T>;



namespace internal
{
  namespace TableBaseAccessors
  {
    template <int N, typename T, bool Constness>
    struct Types
    {};

    template <int N, typename T>
    struct Types<N, T, true>
    {
      using value_type = const T;
      using TableType  = const TableBase<N, T>;

      using iterator       = typename AlignedVector<T>::const_iterator;
      using const_iterator = typename AlignedVector<T>::const_iterator;

      using reference       = typename AlignedVector<T>::const_reference;
      using const_reference = typename AlignedVector<T>::const_reference;
    };

    template <int N, typename T>
    struct Types<N, T, false>
    {
      using value_type = T;
      using TableType  = TableBase<N, T>;

      using iterator       = typename AlignedVector<T>::iterator;
      using const_iterator = typename AlignedVector<T>::const_iterator;

      using reference       = typename AlignedVector<T>::reference;
      using const_reference = typename AlignedVector<T>::const_reference;
    };


    template <int N, typename T, bool C, unsigned int P>
    class Accessor
    {
    public:
      using TableType = typename Types<N, T, C>::TableType;

      using iterator       = typename Types<N, T, C>::iterator;
      using const_iterator = typename Types<N, T, C>::const_iterator;

      using size_type       = size_t;
      using difference_type = ptrdiff_t;

    private:
      Accessor(const TableType &table, const iterator data);

    public:
      Accessor(const Accessor &a);

      Accessor<N, T, C, P - 1> operator[](const size_type i) const;

      DeclException3(ExcIndexRange,
                     size_type,
                     size_type,
                     size_type,
                     << "Index " << N - P + 1 << "has a value of " << arg1
                     << " but needs to be in the range [" << arg2 << "," << arg3
                     << "[.");

    private:
      const TableType &table;
      const iterator   data;

      // declare some other classes
      // as friends. make sure to
      // work around bugs in some
      // compilers
      template <int N1, typename T1>
      friend class ::Table;
      template <int N1, typename T1, bool C1, unsigned int P1>
      friend class Accessor;
#ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
      friend class ::Table<N, T>;
      friend class Accessor<N, T, C, P + 1>;
#endif
    };



    template <int N, typename T, bool C>
    class Accessor<N, T, C, 1>
    {
    public:
      using value_type = typename Types<N, T, C>::value_type;

      using iterator       = typename Types<N, T, C>::iterator;
      using const_iterator = typename Types<N, T, C>::const_iterator;

      using reference       = typename Types<N, T, C>::reference;
      using const_reference = typename Types<N, T, C>::const_reference;

      using size_type       = size_t;
      using difference_type = ptrdiff_t;

      using TableType = typename Types<N, T, C>::TableType;

    private:
      Accessor(const TableType &table, const iterator data);

    public:
      Accessor(const Accessor &a);

      reference operator[](const size_type) const;

      size_type
      size() const;

      iterator
      begin() const;

      iterator
      end() const;

    private:
      const TableType &table;
      const iterator   data;

      // declare some other classes
      // as friends. make sure to
      // work around bugs in some
      // compilers
      template <int N1, typename T1>
      friend class ::Table;
      template <int N1, typename T1, bool C1, unsigned int P1>
      friend class Accessor;
#ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
      friend class ::Table<2, T>;
      friend class Accessor<N, T, C, 2>;
#endif
    };
  } // namespace TableBaseAccessors

} // namespace internal



template <int N, typename T>
class TableBase : public Subscriptor
{
public:
  using value_type = T;

  using size_type = typename AlignedVector<T>::size_type;


  TableBase() = default;

  TableBase(const TableIndices<N> &sizes);

  template <typename InputIterator>
  TableBase(const TableIndices<N> &sizes,
            InputIterator          entries,
            const bool             C_style_indexing = true);

  TableBase(const TableBase<N, T> &src);

  template <typename T2>
  TableBase(const TableBase<N, T2> &src);

  TableBase(TableBase<N, T> &&src) noexcept;

  ~TableBase() override = default;

  TableBase<N, T> &
  operator=(const TableBase<N, T> &src);

  template <typename T2>
  TableBase<N, T> &
  operator=(const TableBase<N, T2> &src);

  TableBase<N, T> &
  operator=(TableBase<N, T> &&src) noexcept;

  bool
  operator==(const TableBase<N, T> &T2) const;

  void
  reset_values();

  void
  reinit(const TableIndices<N> &new_size,
         const bool             omit_default_initialization = false);

  size_type
  size(const unsigned int i) const;

  const TableIndices<N> &
  size() const;

  size_type
  n_elements() const;

  bool
  empty() const;

  template <typename InputIterator>
  void
  fill(InputIterator entries, const bool C_style_indexing = true);

  void
  fill(const T &value);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<N> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<N> &indices) const;

  void
  swap(TableBase<N, T> &v);

  std::size_t
  memory_consumption() const;

  template <class Archive>
  void
  serialize(Archive &ar, const unsigned int version);

protected:
  size_type
  position(const TableIndices<N> &indices) const;

  typename AlignedVector<T>::reference
  el(const TableIndices<N> &indices);

  typename AlignedVector<T>::const_reference
  el(const TableIndices<N> &indices) const;

protected:
  AlignedVector<T> values;

  TableIndices<N> table_size;

  template <int, typename>
  friend class TableBase;
};


template <int N, typename T>
class Table : public TableBase<N, T>
{};


template <typename T>
class Table<1, T> : public TableBase<1, T>
{
public:
  using size_type = typename TableBase<1, T>::size_type;

  Table() = default;

  Table(const size_type size);

  template <typename InputIterator>
  Table(const size_type size,
        InputIterator   entries,
        const bool      C_style_indexing = true);

  typename AlignedVector<T>::const_reference
  operator[](const size_type i) const;

  typename AlignedVector<T>::reference operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i) const;

  typename AlignedVector<T>::reference
  operator()(const size_type i);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<1> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<1> &indices) const;
};



namespace MatrixTableIterators
{
  enum class Storage
  {
    row_major,

    column_major
  };

  // Forward declaration of the iterator class.
  template <typename TableType, bool Constness, Storage storage_order>
  class Iterator;

  template <typename TableType, bool Constness, Storage storage_order>
  class Accessor;

  template <typename TableType, bool Constness, Storage storage_order>
  class AccessorBase
  {
  public:
    using container_pointer_type = typename std::
      conditional<Constness, const TableType *, TableType *>::type;

    using value_type = typename TableType::value_type;

    using size_type = typename TableType::size_type;

    AccessorBase();

    AccessorBase(const container_pointer_type table);

    AccessorBase(const AccessorBase<TableType, false, storage_order> &);

    AccessorBase(const container_pointer_type table,
                 const std::ptrdiff_t         linear_index);

    const value_type &
    value() const;

    size_type
    row() const;

    size_type
    column() const;

  protected:
    container_pointer_type container;

    std::ptrdiff_t linear_index;

    void
    assert_valid_linear_index() const;

    friend class AccessorBase<TableType, true, storage_order>;

    friend class LinearIndexIterator<
      Iterator<TableType, Constness, storage_order>,
      Accessor<TableType, Constness, storage_order>>;
  };

  template <typename TableType, Storage storage_order>
  class Accessor<TableType, true, storage_order>
    : public AccessorBase<TableType, true, storage_order>
  {
  public:
    using value_type =
      typename AccessorBase<TableType, true, storage_order>::value_type;

    using size_type =
      typename AccessorBase<TableType, true, storage_order>::size_type;

    using AccessorBase<TableType, true, storage_order>::AccessorBase;
  };

  template <typename TableType, Storage storage_order>
  class Accessor<TableType, false, storage_order>
    : public AccessorBase<TableType, false, storage_order>
  {
  public:
    using value_type =
      typename AccessorBase<TableType, true, storage_order>::value_type;

    using size_type =
      typename AccessorBase<TableType, true, storage_order>::size_type;

    using AccessorBase<TableType, false, storage_order>::AccessorBase;

    const Accessor<TableType, false, storage_order> &
    operator=(const value_type &) const;

    const Accessor<TableType, false, storage_order> &
    operator=(value_type &&) const;

    using AccessorBase<TableType, false, storage_order>::value;

    value_type &
    value() const;
  };

  template <typename TableType, bool Constness, Storage storage_order>
  class Iterator
    : public LinearIndexIterator<Iterator<TableType, Constness, storage_order>,
                                 Accessor<TableType, Constness, storage_order>>
  {
  public:
    using size_type = typename TableType::size_type;

    using container_pointer_type = typename std::
      conditional<Constness, const TableType *, TableType *>::type;

    Iterator(const Accessor<TableType, Constness, storage_order> &accessor);

    Iterator(const container_pointer_type object);

    Iterator(const container_pointer_type object,
             const size_type              row,
             const size_type              column);

    Iterator(const Iterator<TableType, false, storage_order> &i);

    Iterator(const container_pointer_type container,
             const std::ptrdiff_t         linear_index);
  };
} // namespace MatrixTableIterators



template <typename T>
class Table<2, T> : public TableBase<2, T>
{
public:
  using size_type = typename TableBase<2, T>::size_type;

  using value_type = typename AlignedVector<T>::value_type;

  using reference = typename AlignedVector<T>::reference;

  using const_reference = typename AlignedVector<T>::const_reference;

  using const_iterator = MatrixTableIterators::
    Iterator<Table<2, T>, true, MatrixTableIterators::Storage::row_major>;

  using iterator = MatrixTableIterators::
    Iterator<Table<2, T>, false, MatrixTableIterators::Storage::row_major>;

  Table() = default;

  Table(const size_type size1, const size_type size2);

  template <typename InputIterator>
  Table(const size_type size1,
        const size_type size2,
        InputIterator   entries,
        const bool      C_style_indexing = true);

  void
  reinit(const size_type size1,
         const size_type size2,
         const bool      omit_default_initialization = false);

  using TableBase<2, T>::reinit;

  ::internal::TableBaseAccessors::Accessor<2, T, true, 1>
  operator[](const size_type i) const;

  ::internal::TableBaseAccessors::Accessor<2, T, false, 1>
  operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i, const size_type j) const;


  typename AlignedVector<T>::reference
  operator()(const size_type i, const size_type j);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<2> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<2> &indices) const;


  size_type
  n_rows() const;

  size_type
  n_cols() const;

  iterator
  begin();

  const_iterator
  begin() const;

  iterator
  end();

  const_iterator
  end() const;

protected:
  typename AlignedVector<T>::reference
  el(const size_type i, const size_type j);

  typename AlignedVector<T>::const_reference
  el(const size_type i, const size_type j) const;

  friend class MatrixTableIterators::
    AccessorBase<Table<2, T>, true, MatrixTableIterators::Storage::row_major>;
  friend class MatrixTableIterators::
    AccessorBase<Table<2, T>, false, MatrixTableIterators::Storage::row_major>;

  friend class MatrixTableIterators::
    Accessor<Table<2, T>, false, MatrixTableIterators::Storage::row_major>;
};



template <typename T>
class Table<3, T> : public TableBase<3, T>
{
public:
  using size_type = typename TableBase<3, T>::size_type;

  Table() = default;

  Table(const size_type size1, const size_type size2, const size_type size3);

  template <typename InputIterator>
  Table(const size_type size1,
        const size_type size2,
        const size_type size3,
        InputIterator   entries,
        const bool      C_style_indexing = true);

  ::internal::TableBaseAccessors::Accessor<3, T, true, 2>
  operator[](const size_type i) const;

  ::internal::TableBaseAccessors::Accessor<3, T, false, 2>
  operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i, const size_type j, const size_type k) const;


  typename AlignedVector<T>::reference
  operator()(const size_type i, const size_type j, const size_type k);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<3> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<3> &indices) const;
};



template <typename T>
class Table<4, T> : public TableBase<4, T>
{
public:
  using size_type = typename TableBase<4, T>::size_type;

  Table() = default;

  Table(const size_type size1,
        const size_type size2,
        const size_type size3,
        const size_type size4);

  ::internal::TableBaseAccessors::Accessor<4, T, true, 3>
  operator[](const size_type i) const;

  ::internal::TableBaseAccessors::Accessor<4, T, false, 3>
  operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l) const;


  typename AlignedVector<T>::reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<4> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<4> &indices) const;
};



template <typename T>
class Table<5, T> : public TableBase<5, T>
{
public:
  using size_type = typename TableBase<5, T>::size_type;


  Table() = default;

  Table(const size_type size1,
        const size_type size2,
        const size_type size3,
        const size_type size4,
        const size_type size5);

  ::internal::TableBaseAccessors::Accessor<5, T, true, 4>
  operator[](const size_type i) const;

  ::internal::TableBaseAccessors::Accessor<5, T, false, 4>
  operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l,
             const size_type m) const;

  typename AlignedVector<T>::reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l,
             const size_type m);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<5> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<5> &indices) const;
};



template <typename T>
class Table<6, T> : public TableBase<6, T>
{
public:
  using size_type = typename TableBase<6, T>::size_type;

  Table() = default;

  Table(const size_type size1,
        const size_type size2,
        const size_type size3,
        const size_type size4,
        const size_type size5,
        const size_type size6);

  ::internal::TableBaseAccessors::Accessor<6, T, true, 5>
  operator[](const size_type i) const;

  ::internal::TableBaseAccessors::Accessor<6, T, false, 5>
  operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l,
             const size_type m,
             const size_type n) const;

  typename AlignedVector<T>::reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l,
             const size_type m,
             const size_type n);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<6> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<6> &indices) const;
};


template <typename T>
class Table<7, T> : public TableBase<7, T>
{
public:
  using size_type = typename TableBase<7, T>::size_type;

  Table() = default;

  Table(const size_type size1,
        const size_type size2,
        const size_type size3,
        const size_type size4,
        const size_type size5,
        const size_type size6,
        const size_type size7);

  ::internal::TableBaseAccessors::Accessor<7, T, true, 6>
  operator[](const size_type i) const;

  ::internal::TableBaseAccessors::Accessor<7, T, false, 6>
  operator[](const size_type i);

  typename AlignedVector<T>::const_reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l,
             const size_type m,
             const size_type n,
             const size_type o) const;

  typename AlignedVector<T>::reference
  operator()(const size_type i,
             const size_type j,
             const size_type k,
             const size_type l,
             const size_type m,
             const size_type n,
             const size_type o);

  typename AlignedVector<T>::reference
  operator()(const TableIndices<7> &indices);

  typename AlignedVector<T>::const_reference
  operator()(const TableIndices<7> &indices) const;
};


namespace TransposeTableIterators
{
  template <typename T, bool Constness>
  using AccessorBase DEAL_II_DEPRECATED = MatrixTableIterators::AccessorBase<
    TransposeTable<T>,
    Constness,
    MatrixTableIterators::Storage::column_major>;

  template <typename T, bool Constness>
  using Accessor DEAL_II_DEPRECATED =
    MatrixTableIterators::Accessor<TransposeTable<T>,
                                   Constness,
                                   MatrixTableIterators::Storage::column_major>;

  template <typename T, bool Constness>
  using Iterator DEAL_II_DEPRECATED =
    MatrixTableIterators::Iterator<TransposeTable<T>,
                                   Constness,
                                   MatrixTableIterators::Storage::column_major>;
} // namespace TransposeTableIterators


template <typename T>
class TransposeTable : public TableBase<2, T>
{
public:
  using size_type = typename TableBase<2, T>::size_type;

  using value_type = typename AlignedVector<T>::value_type;

  using reference = typename AlignedVector<T>::reference;

  using const_reference = typename AlignedVector<T>::const_reference;

  using const_iterator =
    MatrixTableIterators::Iterator<TransposeTable<T>,
                                   true,
                                   MatrixTableIterators::Storage::column_major>;

  using iterator =
    MatrixTableIterators::Iterator<TransposeTable<T>,
                                   false,
                                   MatrixTableIterators::Storage::column_major>;

  TransposeTable() = default;

  TransposeTable(const size_type size1, const size_type size2);

  void
  reinit(const size_type size1,
         const size_type size2,
         const bool      omit_default_initialization = false);

  const_reference
  operator()(const size_type i, const size_type j) const;

  reference
  operator()(const size_type i, const size_type j);

  size_type
  n_rows() const;

  size_type
  n_cols() const;

  iterator
  begin();

  const_iterator
  begin() const;

  iterator
  end();

  const_iterator
  end() const;

protected:
  reference
  el(const size_type i, const size_type j);

  const_reference
  el(const size_type i, const size_type j) const;

  friend class MatrixTableIterators::AccessorBase<
    TransposeTable<T>,
    true,
    MatrixTableIterators::Storage::column_major>;
  friend class MatrixTableIterators::AccessorBase<
    TransposeTable<T>,
    false,
    MatrixTableIterators::Storage::column_major>;

  friend class MatrixTableIterators::Accessor<
    TransposeTable<T>,
    false,
    MatrixTableIterators::Storage::column_major>;
};



/* --------------------- Template and inline functions ---------------- */

#ifndef DOXYGEN

template <int N, typename T>
TableBase<N, T>::TableBase(const TableIndices<N> &sizes)
{
  reinit(sizes);
}



template <int N, typename T>
template <typename InputIterator>
TableBase<N, T>::TableBase(const TableIndices<N> &sizes,
                           InputIterator          entries,
                           const bool             C_style_indexing)
{
  reinit(sizes);
  fill(entries, C_style_indexing);
}



template <int N, typename T>
TableBase<N, T>::TableBase(const TableBase<N, T> &src)
  : Subscriptor()
{
  reinit(src.table_size, true);
  values = src.values;
}



template <int N, typename T>
template <typename T2>
TableBase<N, T>::TableBase(const TableBase<N, T2> &src)
{
  reinit(src.table_size);
  if (src.n_elements() != 0)
    std::copy(src.values.begin(), src.values.end(), values.begin());
}



template <int N, typename T>
TableBase<N, T>::TableBase(TableBase<N, T> &&src) noexcept
  : Subscriptor(std::move(src))
  , values(std::move(src.values))
  , table_size(src.table_size)
{
  src.table_size = TableIndices<N>();
}



template <int N, typename T>
template <class Archive>
inline void
TableBase<N, T>::serialize(Archive &ar, const unsigned int)
{
  ar &static_cast<Subscriptor &>(*this);

  ar &values &table_size;
}



namespace internal
{
  namespace TableBaseAccessors
  {
    template <int N, typename T, bool C, unsigned int P>
    inline Accessor<N, T, C, P>::Accessor(const TableType &table,
                                          const iterator   data)
      : table(table)
      , data(data)
    {}



    template <int N, typename T, bool C, unsigned int P>
    inline Accessor<N, T, C, P>::Accessor(const Accessor &a)
      : table(a.table)
      , data(a.data)
    {}



    template <int N, typename T, bool C, unsigned int P>
    inline Accessor<N, T, C, P - 1> Accessor<N, T, C, P>::
                                    operator[](const size_type i) const
    {
      Assert(i < table.size()[N - P], ExcIndexRange(i, 0, table.size()[N - P]));

      // access i-th
      // subobject. optimize on the
      // case i==0
      if (i == 0)
        return Accessor<N, T, C, P - 1>(table, data);
      else
        {
          // note: P>1, otherwise the
          // specialization would have
          // been taken!
          size_type subobject_size = table.size()[N - 1];
          for (int p = P - 1; p > 1; --p)
            subobject_size *= table.size()[N - p];
          const iterator new_data = data + i * subobject_size;
          return Accessor<N, T, C, P - 1>(table, new_data);
        }
    }



    template <int N, typename T, bool C>
    inline Accessor<N, T, C, 1>::Accessor(const TableType &table,
                                          const iterator   data)
      : table(table)
      , data(data)
    {}



    template <int N, typename T, bool C>
    inline Accessor<N, T, C, 1>::Accessor(const Accessor &a)
      : table(a.table)
      , data(a.data)
    {}



    template <int N, typename T, bool C>
    inline typename Accessor<N, T, C, 1>::reference Accessor<N, T, C, 1>::
                                                    operator[](const size_type i) const
    {
      AssertIndexRange(i, table.size()[N - 1]);
      return *(data + i);
    }



    template <int N, typename T, bool C>
    inline typename Accessor<N, T, C, 1>::size_type
    Accessor<N, T, C, 1>::size() const
    {
      return table.size()[N - 1];
    }



    template <int N, typename T, bool C>
    inline typename Accessor<N, T, C, 1>::iterator
    Accessor<N, T, C, 1>::begin() const
    {
      return data;
    }



    template <int N, typename T, bool C>
    inline typename Accessor<N, T, C, 1>::iterator
    Accessor<N, T, C, 1>::end() const
    {
      return data + table.size()[N - 1];
    }
  } // namespace TableBaseAccessors
} // namespace internal



template <int N, typename T>
inline TableBase<N, T> &
TableBase<N, T>::operator=(const TableBase<N, T> &m)
{
  if (!m.empty())
    values = m.values;
  reinit(m.size(), true);

  return *this;
}



template <int N, typename T>
template <typename T2>
inline TableBase<N, T> &
TableBase<N, T>::operator=(const TableBase<N, T2> &m)
{
  reinit(m.size(), true);
  if (!empty())
    std::copy(m.values.begin(),
              m.values.begin() + n_elements(),
              values.begin());

  return *this;
}



template <int N, typename T>
inline TableBase<N, T> &
TableBase<N, T>::operator=(TableBase<N, T> &&m) noexcept
{
  static_cast<Subscriptor &>(*this) = std::move(m);
  values                            = std::move(m.values);
  table_size                        = m.table_size;
  m.table_size                      = TableIndices<N>();

  return *this;
}



template <int N, typename T>
inline bool
TableBase<N, T>::operator==(const TableBase<N, T> &T2) const
{
  return (values == T2.values);
}



template <int N, typename T>
inline void
TableBase<N, T>::reset_values()
{
  // use parallel set operation
  if (n_elements() != 0)
    values.fill(T());
}



template <int N, typename T>
inline void
TableBase<N, T>::fill(const T &value)
{
  if (n_elements() != 0)
    values.fill(value);
}



template <int N, typename T>
inline void
TableBase<N, T>::reinit(const TableIndices<N> &new_sizes,
                        const bool             omit_default_initialization)
{
  table_size = new_sizes;

  const size_type new_size = n_elements();

  // if zero size was given: free all memory
  if (new_size == 0)
    {
      values.resize(0);
      // set all sizes to zero, even
      // if one was previously
      // nonzero. This simplifies
      // some assertions.
      table_size = TableIndices<N>();

      return;
    }

  // adjust values field. If it was empty before, we can simply call resize(),
  // which can set all the data fields. Otherwise, select the fast resize and
  // manually fill in all the elements as required by the design of this
  // class. (Selecting another code for the empty case ensures that we touch
  // the memory only once for non-trivial classes that need to initialize the
  // memory also in resize_fast.)
  if (!omit_default_initialization)
    {
      if (values.empty())
        values.resize(new_size);
      else
        {
          values.resize_fast(new_size);
          values.fill();
        }
    }
  else
    values.resize_fast(new_size);
}



template <int N, typename T>
inline const TableIndices<N> &
TableBase<N, T>::size() const
{
  return table_size;
}



template <int N, typename T>
inline typename TableBase<N, T>::size_type
TableBase<N, T>::size(const unsigned int i) const
{
  Assert(i < N, ExcIndexRange(i, 0, N));
  return table_size[i];
}



template <int N, typename T>
inline typename TableBase<N, T>::size_type
TableBase<N, T>::n_elements() const
{
  size_type s = 1;
  for (unsigned int n = 0; n < N; ++n)
    s *= table_size[n];
  return s;
}



template <int N, typename T>
inline bool
TableBase<N, T>::empty() const
{
  return (n_elements() == 0);
}



namespace internal
{
  namespace TableImplementation
  {
    template <typename InputIterator, typename T>
    void
    fill_Fortran_style(InputIterator entries, TableBase<1, T> &table)
    {
      using size_type = typename TableBase<1, T>::size_type;
      for (size_type i = 0; i < table.size()[0]; ++i)
        table(TableIndices<1>(i)) = *entries++;
    }


    template <typename InputIterator, typename T>
    void
    fill_Fortran_style(InputIterator entries, TableBase<2, T> &table)
    {
      using size_type = typename TableBase<2, T>::size_type;
      for (size_type j = 0; j < table.size()[1]; ++j)
        for (size_type i = 0; i < table.size()[0]; ++i)
          table(TableIndices<2>(i, j)) = *entries++;
    }


    template <typename InputIterator, typename T>
    void
    fill_Fortran_style(InputIterator entries, TableBase<3, T> &table)
    {
      using size_type = typename TableBase<3, T>::size_type;
      for (size_type k = 0; k < table.size()[2]; ++k)
        for (size_type j = 0; j < table.size()[1]; ++j)
          for (size_type i = 0; i < table.size()[0]; ++i)
            table(TableIndices<3>(i, j, k)) = *entries++;
    }


    template <typename InputIterator, typename T, int N>
    void
    fill_Fortran_style(InputIterator, TableBase<N, T> &)
    {
      Assert(false, ExcNotImplemented());
    }
  } // namespace TableImplementation
} // namespace internal


template <int N, typename T>
template <typename InputIterator>
inline void
TableBase<N, T>::fill(InputIterator entries, const bool C_style_indexing)
{
  Assert(n_elements() != 0, ExcMessage("Trying to fill an empty matrix."));

  if (C_style_indexing)
    for (typename AlignedVector<T>::iterator p = values.begin();
         p != values.end();
         ++p)
      *p = *entries++;
  else
    internal::TableImplementation::fill_Fortran_style(entries, *this);
}



template <int N, typename T>
inline void
TableBase<N, T>::swap(TableBase<N, T> &v)
{
  values.swap(v.values);
  std::swap(table_size, v.table_size);
}



template <int N, typename T>
inline std::size_t
TableBase<N, T>::memory_consumption() const
{
  return sizeof(*this) + MemoryConsumption::memory_consumption(values);
}



template <int N, typename T>
inline typename TableBase<N, T>::size_type
TableBase<N, T>::position(const TableIndices<N> &indices) const
{
  // specialize this for the
  // different numbers of dimensions,
  // to make the job somewhat easier
  // for the compiler. have the
  // general formula nevertheless:
  switch (N)
    {
      case 1:
        return indices[0];
      case 2:
        return size_type(indices[0]) * table_size[1] + indices[1];
      case 3:
        return ((size_type(indices[0]) * table_size[1] + indices[1]) *
                  table_size[2] +
                indices[2]);
      default:
        {
          unsigned int s = indices[0];
          for (unsigned int n = 1; n < N; ++n)
            s = s * table_size[n] + indices[n];
          return s;
        }
    }
}



template <int N, typename T>
inline typename AlignedVector<T>::const_reference
TableBase<N, T>::operator()(const TableIndices<N> &indices) const
{
  for (unsigned int n = 0; n < N; ++n)
    AssertIndexRange(indices[n], table_size[n]);
  return el(indices);
}



template <int N, typename T>
inline typename AlignedVector<T>::reference
TableBase<N, T>::operator()(const TableIndices<N> &indices)
{
  for (unsigned int n = 0; n < N; ++n)
    AssertIndexRange(indices[n], table_size[n]);
  return el(indices);
}



template <int N, typename T>
inline typename AlignedVector<T>::const_reference
TableBase<N, T>::el(const TableIndices<N> &indices) const
{
  return values[position(indices)];
}



template <int N, typename T>
inline typename AlignedVector<T>::reference
TableBase<N, T>::el(const TableIndices<N> &indices)
{
  AssertIndexRange(position(indices), values.size());
  return values[position(indices)];
}



template <typename T>
inline Table<1, T>::Table(const size_type size)
  : TableBase<1, T>(TableIndices<1>(size))
{}



template <typename T>
template <typename InputIterator>
inline Table<1, T>::Table(const size_type size,
                          InputIterator   entries,
                          const bool      C_style_indexing)
  : TableBase<1, T>(TableIndices<1>(size), entries, C_style_indexing)
{}



template <typename T>
inline typename AlignedVector<T>::const_reference Table<1, T>::
                                                  operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  return this->values[i];
}



template <typename T>
inline typename AlignedVector<T>::reference Table<1, T>::
                                            operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  return this->values[i];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<1, T>::operator()(const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  return this->values[i];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<1, T>::operator()(const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  return this->values[i];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<1, T>::operator()(const TableIndices<1> &indices) const
{
  return TableBase<1, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<1, T>::operator()(const TableIndices<1> &indices)
{
  return TableBase<1, T>::operator()(indices);
}


//---------------------------------------------------------------------------



template <typename T>
inline Table<2, T>::Table(const size_type size1, const size_type size2)
  : TableBase<2, T>(TableIndices<2>(size1, size2))
{}



template <typename T>
template <typename InputIterator>
inline Table<2, T>::Table(const size_type size1,
                          const size_type size2,
                          InputIterator   entries,
                          const bool      C_style_indexing)
  : TableBase<2, T>(TableIndices<2>(size1, size2), entries, C_style_indexing)
{}



template <typename T>
inline void
Table<2, T>::reinit(const size_type size1,
                    const size_type size2,
                    const bool      omit_default_initialization)
{
  this->TableBase<2, T>::reinit(TableIndices<2>(size1, size2),
                                omit_default_initialization);
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<2, T, true, 1>
  Table<2, T>::operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  return ::internal::TableBaseAccessors::Accessor<2, T, true, 1>(
    *this, this->values.begin() + size_type(i) * n_cols());
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<2, T, false, 1>
  Table<2, T>::operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  return ::internal::TableBaseAccessors::Accessor<2, T, false, 1>(
    *this, this->values.begin() + size_type(i) * n_cols());
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<2, T>::operator()(const size_type i, const size_type j) const
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  return this->values[size_type(i) * this->table_size[1] + j];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<2, T>::operator()(const size_type i, const size_type j)
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  return this->values[size_type(i) * this->table_size[1] + j];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<2, T>::operator()(const TableIndices<2> &indices) const
{
  return TableBase<2, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<2, T>::operator()(const TableIndices<2> &indices)
{
  return TableBase<2, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<2, T>::el(const size_type i, const size_type j) const
{
  return this->values[size_type(i) * this->table_size[1] + j];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<2, T>::el(const size_type i, const size_type j)
{
  return this->values[size_type(i) * this->table_size[1] + j];
}



template <typename T>
inline typename Table<2, T>::size_type
Table<2, T>::n_rows() const
{
  return this->table_size[0];
}



template <typename T>
inline typename Table<2, T>::size_type
Table<2, T>::n_cols() const
{
  return this->table_size[1];
}



template <typename T>
inline typename Table<2, T>::iterator
Table<2, T>::begin()
{
  return typename Table<2, T>::iterator(this, 0, 0);
}



template <typename T>
inline typename Table<2, T>::const_iterator
Table<2, T>::begin() const
{
  return typename Table<2, T>::const_iterator(this, 0, 0);
}



template <typename T>
inline typename Table<2, T>::iterator
Table<2, T>::end()
{
  return typename Table<2, T>::iterator(this);
}



template <typename T>
inline typename Table<2, T>::const_iterator
Table<2, T>::end() const
{
  return typename Table<2, T>::const_iterator(this);
}



//---------------------------------------------------------------------------
namespace MatrixTableIterators
{
  namespace internal
  {
    // Internal calculation routines for AccessorBase. These do not do any
    // checking whatsoever.
    template <typename TableType, Storage storage_order>
    inline std::ptrdiff_t
    get_row_index(const std::ptrdiff_t   linear_index,
                  const TableType *const container)
    {
      switch (storage_order)
        {
          case Storage::row_major:
            return linear_index / container->n_cols();
          case Storage::column_major:
            return linear_index % container->n_rows();
          default:
            Assert(false, ExcInternalError());
        }
      return {};
    }



    template <typename TableType, Storage storage_order>
    inline std::ptrdiff_t
    get_column_index(const std::ptrdiff_t   linear_index,
                     const TableType *const container)
    {
      switch (storage_order)
        {
          case Storage::row_major:
            return linear_index % container->n_cols();
          case Storage::column_major:
            return linear_index / container->n_rows();
          default:
            Assert(false, ExcInternalError());
        }
      return {};
    }
  } // namespace internal



  template <typename TableType, bool Constness, Storage storage_order>
  inline AccessorBase<TableType, Constness, storage_order>::AccessorBase()
    : container(nullptr)
    , linear_index(std::numeric_limits<decltype(linear_index)>::max())
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline AccessorBase<TableType, Constness, storage_order>::AccessorBase(
    const container_pointer_type table)
    : container(table)
    , linear_index(container->values.size())
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline AccessorBase<TableType, Constness, storage_order>::AccessorBase(
    const AccessorBase<TableType, false, storage_order> &a)
    : container(a.container)
    , linear_index(a.linear_index)
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline AccessorBase<TableType, Constness, storage_order>::AccessorBase(
    const container_pointer_type table,
    const std::ptrdiff_t         index)
    : container(table)
    , linear_index(index)
  {
    Assert(0 <= linear_index &&
             std::size_t(linear_index) < container->values.size() + 1,
           ExcMessage("The current iterator points outside of the table and is "
                      "not the end iterator."));
  }



  template <typename TableType, bool Constness, Storage storage_order>
  inline const typename AccessorBase<TableType, Constness, storage_order>::
    value_type &
    AccessorBase<TableType, Constness, storage_order>::value() const
  {
    assert_valid_linear_index();
    return this->container->values[linear_index];
  }



  template <typename TableType, bool Constness, Storage storage_order>
  inline typename AccessorBase<TableType, Constness, storage_order>::size_type
  AccessorBase<TableType, Constness, storage_order>::row() const
  {
    assert_valid_linear_index();
    return static_cast<std::size_t>(
      internal::get_row_index<TableType, storage_order>(linear_index,
                                                        container));
  }



  template <typename TableType, bool Constness, Storage storage_order>
  inline typename AccessorBase<TableType, Constness, storage_order>::size_type
  AccessorBase<TableType, Constness, storage_order>::column() const
  {
    assert_valid_linear_index();
    return static_cast<std::size_t>(
      internal::get_column_index<TableType, storage_order>(linear_index,
                                                           container));
  }



  template <typename TableType, bool Constness, Storage storage_order>
  inline void
  AccessorBase<TableType, Constness, storage_order>::assert_valid_linear_index()
    const
  {
#  ifdef DEBUG // avoid unused variable warnings by guarding everything
    Assert(container != nullptr,
           ExcMessage("This accessor has been default-constructed and does not "
                      "have a corresponding table."));
    Assert(!container->empty(),
           ExcMessage("An empty table has no rows or columns."));
    Assert(0 <= linear_index &&
             std::size_t(linear_index) < container->values.size(),
           ExcMessage("The current iterator points outside of the table."));
    const std::ptrdiff_t row_n =
      internal::get_row_index<TableType, storage_order>(linear_index,
                                                        container);
    const std::ptrdiff_t column_n =
      internal::get_column_index<TableType, storage_order>(linear_index,
                                                           container);
    Assert(0 <= column_n && std::size_t(column_n) < container->n_cols(),
           ExcMessage("The current iterator points outside the table."));
    Assert(0 <= row_n && std::size_t(row_n) < container->n_rows(),
           ExcMessage("The current iterator points outside the table."));
#  endif
  }



  template <typename TableType, Storage storage_order>
  inline const Accessor<TableType, false, storage_order> &
  Accessor<TableType, false, storage_order>::operator=(
    const typename Accessor<TableType, false, storage_order>::value_type &t)
    const
  {
    this->assert_valid_linear_index();
    this->container->values[this->linear_index] = t;
    return *this;
  }



  template <typename TableType, Storage storage_order>
  inline const Accessor<TableType, false, storage_order> &
  Accessor<TableType, false, storage_order>::operator=(
    typename Accessor<TableType, false, storage_order>::value_type &&t) const
  {
    this->assert_valid_linear_index();
    this->container->values[this->linear_index] = t;
    return *this;
  }



  template <typename TableType, Storage storage_order>
  inline typename Accessor<TableType, false, storage_order>::value_type &
  Accessor<TableType, false, storage_order>::value() const
  {
    this->assert_valid_linear_index();
    return this->container->values[this->linear_index];
  }



  template <typename TableType, bool Constness, Storage storage_order>
  inline Iterator<TableType, Constness, storage_order>::Iterator(
    const Accessor<TableType, Constness, storage_order> &a)
    : LinearIndexIterator<Iterator<TableType, Constness, storage_order>,
                          Accessor<TableType, Constness, storage_order>>(a)
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline Iterator<TableType, Constness, storage_order>::Iterator(
    const container_pointer_type table)
    : LinearIndexIterator<Iterator<TableType, Constness, storage_order>,
                          Accessor<TableType, Constness, storage_order>>(
        Accessor<TableType, Constness, storage_order>(table))
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline Iterator<TableType, Constness, storage_order>::Iterator(
    const Iterator<TableType, false, storage_order> &i)
    : LinearIndexIterator<Iterator<TableType, Constness, storage_order>,
                          Accessor<TableType, Constness, storage_order>>(*i)
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline Iterator<TableType, Constness, storage_order>::Iterator(
    const container_pointer_type table,
    const size_type              row_n,
    const size_type              col_n)
    : Iterator(table,
               storage_order == Storage::row_major ?
                 table->n_cols() * row_n + col_n :
                 table->n_rows() * col_n + row_n)
  {}



  template <typename TableType, bool Constness, Storage storage_order>
  inline Iterator<TableType, Constness, storage_order>::Iterator(
    const container_pointer_type table,
    const std::ptrdiff_t         linear_index)
    : LinearIndexIterator<Iterator<TableType, Constness, storage_order>,
                          Accessor<TableType, Constness, storage_order>>(
        Accessor<TableType, Constness, storage_order>(table, linear_index))
  {}
} // namespace MatrixTableIterators



//---------------------------------------------------------------------------
template <typename T>
inline TransposeTable<T>::TransposeTable(const size_type size1,
                                         const size_type size2)
  : TableBase<2, T>(TableIndices<2>(size2, size1))
{}



template <typename T>
inline void
TransposeTable<T>::reinit(const size_type size1,
                          const size_type size2,
                          const bool      omit_default_initialization)
{
  this->TableBase<2, T>::reinit(TableIndices<2>(size2, size1),
                                omit_default_initialization);
}



template <typename T>
inline typename TransposeTable<T>::const_reference
TransposeTable<T>::operator()(const size_type i, const size_type j) const
{
  AssertIndexRange(i, this->table_size[1]);
  AssertIndexRange(j, this->table_size[0]);
  return this->values[size_type(j) * this->table_size[1] + i];
}



template <typename T>
inline typename TransposeTable<T>::reference
TransposeTable<T>::operator()(const size_type i, const size_type j)
{
  AssertIndexRange(i, this->table_size[1]);
  AssertIndexRange(j, this->table_size[0]);
  return this->values[size_type(j) * this->table_size[1] + i];
}



template <typename T>
inline typename TransposeTable<T>::const_reference
TransposeTable<T>::el(const size_type i, const size_type j) const
{
  return this->values[size_type(j) * this->table_size[1] + i];
}



template <typename T>
inline typename TransposeTable<T>::reference
TransposeTable<T>::el(const size_type i, const size_type j)
{
  return this->values[size_type(j) * this->table_size[1] + i];
}



template <typename T>
inline typename TransposeTable<T>::size_type
TransposeTable<T>::n_rows() const
{
  return this->table_size[1];
}



template <typename T>
inline typename TransposeTable<T>::size_type
TransposeTable<T>::n_cols() const
{
  return this->table_size[0];
}



template <typename T>
inline typename TransposeTable<T>::iterator
TransposeTable<T>::begin()
{
  return typename TransposeTable<T>::iterator(this, 0, 0);
}



template <typename T>
inline typename TransposeTable<T>::const_iterator
TransposeTable<T>::begin() const
{
  return typename TransposeTable<T>::const_iterator(this, 0, 0);
}



template <typename T>
inline typename TransposeTable<T>::iterator
TransposeTable<T>::end()
{
  return typename TransposeTable<T>::iterator(this);
}



template <typename T>
inline typename TransposeTable<T>::const_iterator
TransposeTable<T>::end() const
{
  return typename TransposeTable<T>::const_iterator(this);
}
//---------------------------------------------------------------------------



template <typename T>
inline Table<3, T>::Table(const size_type size1,
                          const size_type size2,
                          const size_type size3)
  : TableBase<3, T>(TableIndices<3>(size1, size2, size3))
{}



template <typename T>
template <typename InputIterator>
inline Table<3, T>::Table(const size_type size1,
                          const size_type size2,
                          const size_type size3,
                          InputIterator   entries,
                          const bool      C_style_indexing)
  : TableBase<3, T>(TableIndices<3>(size1, size2, size3),
                    entries,
                    C_style_indexing)
{}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<3, T, true, 2>
  Table<3, T>::operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size =
    size_type(this->table_size[1]) * this->table_size[2];
  return (::internal::TableBaseAccessors::Accessor<3, T, true, 2>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<3, T, false, 2>
  Table<3, T>::operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size =
    size_type(this->table_size[1]) * this->table_size[2];
  return (::internal::TableBaseAccessors::Accessor<3, T, false, 2>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<3, T>::
operator()(const size_type i, const size_type j, const size_type k) const
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  return this
    ->values[(size_type(i) * this->table_size[1] + j) * this->table_size[2] +
             k];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<3, T>::operator()(const size_type i, const size_type j, const size_type k)
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  return this
    ->values[(size_type(i) * this->table_size[1] + j) * this->table_size[2] +
             k];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<3, T>::operator()(const TableIndices<3> &indices) const
{
  return TableBase<3, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<3, T>::operator()(const TableIndices<3> &indices)
{
  return TableBase<3, T>::operator()(indices);
}



template <typename T>
inline Table<4, T>::Table(const size_type size1,
                          const size_type size2,
                          const size_type size3,
                          const size_type size4)
  : TableBase<4, T>(TableIndices<4>(size1, size2, size3, size4))
{}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<4, T, true, 3>
  Table<4, T>::operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size =
    size_type(this->table_size[1]) * this->table_size[2] * this->table_size[3];
  return (::internal::TableBaseAccessors::Accessor<4, T, true, 3>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<4, T, false, 3>
  Table<4, T>::operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size =
    size_type(this->table_size[1]) * this->table_size[2] * this->table_size[3];
  return (::internal::TableBaseAccessors::Accessor<4, T, false, 3>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<4, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l) const
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  return this
    ->values[((size_type(i) * this->table_size[1] + j) * this->table_size[2] +
              k) *
               this->table_size[3] +
             l];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<4, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l)
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  return this
    ->values[((size_type(i) * this->table_size[1] + j) * this->table_size[2] +
              k) *
               this->table_size[3] +
             l];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<4, T>::operator()(const TableIndices<4> &indices) const
{
  return TableBase<4, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<4, T>::operator()(const TableIndices<4> &indices)
{
  return TableBase<4, T>::operator()(indices);
}



template <typename T>
inline Table<5, T>::Table(const size_type size1,
                          const size_type size2,
                          const size_type size3,
                          const size_type size4,
                          const size_type size5)
  : TableBase<5, T>(TableIndices<5>(size1, size2, size3, size4, size5))
{}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<5, T, true, 4>
  Table<5, T>::operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size = size_type(this->table_size[1]) *
                                   this->table_size[2] * this->table_size[3] *
                                   this->table_size[4];
  return (::internal::TableBaseAccessors::Accessor<5, T, true, 4>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<5, T, false, 4>
  Table<5, T>::operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size = size_type(this->table_size[1]) *
                                   this->table_size[2] * this->table_size[3] *
                                   this->table_size[4];
  return (::internal::TableBaseAccessors::Accessor<5, T, false, 4>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<5, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l,
                        const size_type m) const
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  AssertIndexRange(m, this->table_size[4]);
  return this
    ->values[(((size_type(i) * this->table_size[1] + j) * this->table_size[2] +
               k) *
                this->table_size[3] +
              l) *
               this->table_size[4] +
             m];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<5, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l,
                        const size_type m)
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  AssertIndexRange(m, this->table_size[4]);
  return this
    ->values[(((size_type(i) * this->table_size[1] + j) * this->table_size[2] +
               k) *
                this->table_size[3] +
              l) *
               this->table_size[4] +
             m];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<5, T>::operator()(const TableIndices<5> &indices) const
{
  return TableBase<5, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<5, T>::operator()(const TableIndices<5> &indices)
{
  return TableBase<5, T>::operator()(indices);
}



template <typename T>
inline Table<6, T>::Table(const size_type size1,
                          const size_type size2,
                          const size_type size3,
                          const size_type size4,
                          const size_type size5,
                          const size_type size6)
  : TableBase<6, T>()
{
  TableIndices<6> table_indices;
  table_indices[0] = size1;
  table_indices[1] = size2;
  table_indices[2] = size3;
  table_indices[3] = size4;
  table_indices[4] = size5;
  table_indices[5] = size6;

  TableBase<6, T>::reinit(table_indices);
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<6, T, true, 5>
  Table<6, T>::operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size = size_type(this->table_size[1]) *
                                   this->table_size[2] * this->table_size[3] *
                                   this->table_size[4] * this->table_size[5];
  return (::internal::TableBaseAccessors::Accessor<6, T, true, 5>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<6, T, false, 5>
  Table<6, T>::operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size = size_type(this->table_size[1]) *
                                   this->table_size[2] * this->table_size[3] *
                                   this->table_size[4] * this->table_size[5];
  return (::internal::TableBaseAccessors::Accessor<6, T, false, 5>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<6, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l,
                        const size_type m,
                        const size_type n) const
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  AssertIndexRange(m, this->table_size[4]);
  AssertIndexRange(n, this->table_size[5]);
  return this
    ->values[((((size_type(i) * this->table_size[1] + j) * this->table_size[2] +
                k) *
                 this->table_size[3] +
               l) *
                this->table_size[4] +
              m) *
               this->table_size[5] +
             n];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<6, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l,
                        const size_type m,
                        const size_type n)
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  AssertIndexRange(m, this->table_size[4]);
  AssertIndexRange(n, this->table_size[5]);
  return this
    ->values[((((size_type(i) * this->table_size[1] + j) * this->table_size[2] +
                k) *
                 this->table_size[3] +
               l) *
                this->table_size[4] +
              m) *
               this->table_size[5] +
             n];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<6, T>::operator()(const TableIndices<6> &indices) const
{
  return TableBase<6, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<6, T>::operator()(const TableIndices<6> &indices)
{
  return TableBase<6, T>::operator()(indices);
}



template <typename T>
inline Table<7, T>::Table(const size_type size1,
                          const size_type size2,
                          const size_type size3,
                          const size_type size4,
                          const size_type size5,
                          const size_type size6,
                          const size_type size7)
  : TableBase<7, T>()
{
  TableIndices<7> table_indices;
  table_indices[0] = size1;
  table_indices[1] = size2;
  table_indices[2] = size3;
  table_indices[3] = size4;
  table_indices[4] = size5;
  table_indices[5] = size6;
  table_indices[6] = size7;

  TableBase<7, T>::reinit(table_indices);
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<7, T, true, 6>
  Table<7, T>::operator[](const size_type i) const
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size =
    size_type(this->table_size[1]) * this->table_size[2] * this->table_size[3] *
    this->table_size[4] * this->table_size[5] * this->table_size[6];
  return (::internal::TableBaseAccessors::Accessor<7, T, true, 6>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline ::internal::TableBaseAccessors::Accessor<7, T, false, 6>
  Table<7, T>::operator[](const size_type i)
{
  AssertIndexRange(i, this->table_size[0]);
  const size_type subobject_size =
    size_type(this->table_size[1]) * this->table_size[2] * this->table_size[3] *
    this->table_size[4] * this->table_size[5] * this->table_size[6];
  return (::internal::TableBaseAccessors::Accessor<7, T, false, 6>(
    *this, this->values.begin() + i * subobject_size));
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<7, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l,
                        const size_type m,
                        const size_type n,
                        const size_type o) const
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  AssertIndexRange(m, this->table_size[4]);
  AssertIndexRange(n, this->table_size[5]);
  AssertIndexRange(o, this->table_size[6]);
  return this->values
    [(((((size_type(i) * this->table_size[1] + j) * this->table_size[2] + k) *
          this->table_size[3] +
        l) *
         this->table_size[4] +
       m) *
        this->table_size[5] +
      n) *
       this->table_size[6] +
     o];
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<7, T>::operator()(const size_type i,
                        const size_type j,
                        const size_type k,
                        const size_type l,
                        const size_type m,
                        const size_type n,
                        const size_type o)
{
  AssertIndexRange(i, this->table_size[0]);
  AssertIndexRange(j, this->table_size[1]);
  AssertIndexRange(k, this->table_size[2]);
  AssertIndexRange(l, this->table_size[3]);
  AssertIndexRange(m, this->table_size[4]);
  AssertIndexRange(n, this->table_size[5]);
  AssertIndexRange(o, this->table_size[6]);
  return this->values
    [(((((size_type(i) * this->table_size[1] + j) * this->table_size[2] + k) *
          this->table_size[3] +
        l) *
         this->table_size[4] +
       m) *
        this->table_size[5] +
      n) *
       this->table_size[6] +
     o];
}



template <typename T>
inline typename AlignedVector<T>::const_reference
Table<7, T>::operator()(const TableIndices<7> &indices) const
{
  return TableBase<7, T>::operator()(indices);
}



template <typename T>
inline typename AlignedVector<T>::reference
Table<7, T>::operator()(const TableIndices<7> &indices)
{
  return TableBase<7, T>::operator()(indices);
}


#endif // DOXYGEN



template <int N, typename T>
inline void
swap(TableBase<N, T> &u, TableBase<N, T> &v)
{
  u.swap(v);
}

DEAL_II_NAMESPACE_CLOSE

#endif