⛏️

C++で行列クラスを作成する(1)

15 min read

C++で行列クラスを作成する。

要件としてはC++20でconstexpr関数内で使用でき、<memory_resource>ヘッダのstd::pmr::polymorphic_allocatorにもしっかり対応している事。

<memory_resource>はVC++しかまだ対応しているのを見ていないのでVC++でC++のバージョンをC++20にして作成するものとする。

長いので記事は分割します。

matrix_type

行列で使用する型で満たしてほしい要件をmatrix_typeコンセプトで定義する。
とりあえず四則演算と代入ができればOKとする。

#include <concepts>
template<typename T>
concept matrix_type = requires(T& a, T& b){
  { a + b } -> std::convertible_to<T>;
  { a - b } -> std::convertible_to<T>;
  { a * b } -> std::convertible_to<T>;
  { a / b } -> std::convertible_to<T>;
  a = b;
  a = std::move(b);
};
static_assert(matrix_type<int>);
static_assert(matrix_type<float>);
static_assert(matrix_type<double>);

必要に応じて要件を追加していくものとする。

dynamic_matrix

要素数は可変になってほしいのでdynamic_matrixという名前にする。

#include <memory>
#include <cassert>

template<matrix_type T, typename Alloc=std::allocator<T>>
class dynamic_matrix{
public:
  using allocator_type  = std::allocator_traits<Alloc>::template rebind_alloc<T>;
  using value_type      = std::allocator_traits<allocator_type>::value_type;
  using pointer         = std::allocator_traits<allocator_type>::pointer;
  using const_pointer   = std::allocator_traits<allocator_type>::const_pointer;
  using reference       = value_type &;
  using const_reference = const value_type &;
  using size_type       = std::allocator_traits<allocator_type>::size_type;
  using index_type      = size_type;
  //メモリアクセス
  constexpr reference operator()(index_type i, index_type j) &noexcept{
    assert(m_elem != nullptr);
    assert(i < m_rows);
    assert(j < m_cols);
    return m_elem[m_colcap * i + j];
  }
  constexpr const_reference operator()(index_type i, index_type j) const &noexcept{
    assert(m_elem != nullptr);
    assert(i < m_rows);
    assert(j < m_cols);
    return m_elem[m_colcap * i + j];
  }
  constexpr value_type operator()(index_type i, index_type j) const &&noexcept{
    assert(m_elem != nullptr);
    assert(i < m_rows);
    assert(j < m_cols);
    return m_elem[m_colcap * i + j];
  }
  //実際に使うサイズ
  constexpr size_type rows()const noexcept {
    return m_rows;
  }
  constexpr size_type columns()const noexcept {
    return m_cols;
  }
  constexpr size_type size()const noexcept {
    return rows() * columns();
  }
  //余分に確保しているメモリサイズ
  constexpr size_type row_capacity()const noexcept {
    return m_rowcap;
  }
  constexpr size_type column_capacity()const noexcept {
    return m_colcap;
  }
  constexpr size_type capacity()const noexcept {
    return row_capacity() * column_capacity();
  }
  //アロケータを返す
  constexpr allocator_type get_allocator() const noexcept {
    return m_alloc;
  }
private:
  using traits = std::allocator_traits<allocator_type>;
  //m_elemがnullptrであるときallocatorされていない、
  //もしくはdeallocate済みで有ることを保証すること。
  //関数から抜けたときにm_elemがnullptr以外であれば
  //allocateされているという事も保証する事
  //constexpr対応の為に必須の要件とする。
  pointer m_elem         = nullptr;
  allocator_type m_alloc = {};
  size_type m_rows       = 0;
  size_type m_cols       = 0;
  size_type m_rowcap     = 0;
  size_type m_colcap     = 0;
};

dynamic_matrixのメモリは以下のようにする予定である。
dynamic_matrixのメモリのイメージ
figure1-1: dynamic_matrixのメモリのイメージ

上記の灰色の部分は余分に確保しているメモリ(***_capacity())で数字の部分は実際に行列として使用する範囲である。
灰色の部分の値にアクセスは未定義とする。

行のインデクスをi、列のインデックスをjとしたときM_{ij}にアクセスしたい場合以下の式となる。

m_elem[m_colcap * i + j]

数学の行列と異なりインデックスは0スタートとなる事には注意
たとえば上記の行列でi=1, j=26となる。
メモリアクセス
figure1-2: メモリアクセス

メモリ確保

上記内容でメモリを確保する関数を作成していこう。
関数名はstd::vectorを参考にする。

public:
//dynamic_matrixのメンバー関数

//メモリを確保する。現在のメモリより小さくはしない。
constexpr void reserve(
  size_type rows, size_type cols,
  size_type rowcap, size_type colcap){
  //rowcap,colcapはそれぞれrows,colsより小さくてはいけない
  if(rows > rowcap)
    rowcap = rows;
  if(cols > colcap)
    colcap = cols;
  //現在のcapacityより大きくメモリを確保したいときのみメモリ確保を行う
  if (row_capacity() < rowcap || column_capacity() < colcap) {
    //一時的なメモリを確保
    pointer tmp = traits::allocate(m_alloc, rowcap * colcap);
   
    //constexpr関数内で使う為にメモリ全て初期化する。
    for (pointer sptr = tmp, eptr = tmp + rowcap * colcap; sptr < eptr; ++sptr)
      traits::construct(m_alloc, sptr);
 
    //すでにm_elemでメモリを確保されている場合にのみtmpにデータをコピー
    if (m_elem) {
      size_type dst_rows = std::min(rows, m_rows);
      size_type dst_cols = std::min(cols, m_cols);

      for (size_type ri = 0; ri < dst_rows; ++ri)
        for (size_type ci = 0; ci < dst_cols; ++ci)
          tmp[colcap * ri + ci] = (*this)(ri, ci);
      traits::deallocate(m_alloc, m_elem, capacity());
    }
    m_elem = tmp;
    m_rowcap = rowcap;
    m_colcap = colcap;
  }

  m_rows   = rows;
  m_cols   = cols;
}
//容量のみの変更を想定
constexpr void reserve(size_type rowcap, size_type colcap){
  reserve(m_rows, m_cols, rowcap, colcap);
}
//サイズの変更を想定。
//ただし現在の容量を超えるサイズであればメモリを再確保する。
constexpr void resize(size_type rows, size_type cols){
  reserve(rows, cols, m_rowcap, m_colcap);
}
//rows()*columns()の範囲を指定の値で初期化
constexpr void fill(const T& value) {
  for (size_type ri = 0; ri < m_rows; ++ri)
    for (size_type ci = 0; ci < m_cols; ++ci)
      (*this)(ri, ci) = value;
}
//デストラクタもconstexprとする事。
constexpr ~dynamic_matrix(){
  deallocate();
}
private:
constexpr void deallocate(bool post_processing=true){
  //重要: m_elemがnullptrの場合、allocateされていないか、
  //すでにdeallocateされているものとする。
  if (m_elem) {
    //後処理
    for (pointer sptr = m_elem, eptr = m_elem + capacity(); sptr < eptr; ++sptr)
      traits::destroy(m_alloc, sptr);
    //メモリ解放
    traits::deallocate(m_alloc, m_elem, capacity());
    m_elem = nullptr;
  }
  //後始末をしないほうがいいときがあるのでpost_processingで分岐する。
  //基本は後始末をしたほうがいいのでデフォルトではtrueである。
  if(post_processing){
    m_elem   = nullptr;
    m_rows   = 0;
    m_cols   = 0;
    m_rowcap = 0;
    m_colcap = 0;
  }
}

初期化(std::initializer_list)

メモリ確保ができるようになったのでコンストラクタで初期化できるようにしたい。

public:
//dynamic_matrixのコンストラクタ
//std::initializer_list<U>のサイズは全て同じである事を条件とする。
//異なる場合未定義とする
template<std::convertible_to<T> U>
constexpr dynamic_matrix(std::initializer_list<std::initializer_list<U>> mat,
  const allocator_type &alloc = allocator_type()):
  //allocatorに代入演算子が存在しない事もあるので初期化はここでのみ。
  m_alloc(alloc){
  size_type M = mat.size();
  size_type N = mat.begin()->size();
  resize(M, N);
  index_type ri = 0;
  for (auto &&row : mat) {
    index_type ci = 0;
    for (auto &&value : row)
      (*this)(ri,ci++) = static_cast<T>(value);
    ++ri;
  }
}

上記コンストラクタで以下のような初期化が可能となる。

dynamic_matrix<double> mat({
  {1, 2, 3},
  {4, 5, 6},
  {7, 8, 9}
});

出力

初期化ができるようになったなら出力もできるようにしたほうがわかりやすい。

#include <iostream>
#include <iomanip>

template <typename Char, matrix_type T, typename Alloc>
std::basic_ostream<Char> &operator<<(
  std::basic_ostream<Char> &os, const dynamic_matrix<T, Alloc> &mat) {
  using size_type = dynamic_matrix<T, Alloc>::size_type;
  size_type M     = mat.rows();
  size_type N     = mat.columns();
  os << std::showpos << std::fixed << std::setprecision(3);
  for (size_type ri = 0; ri < M; ++ri) {
    for (size_type ci = 0; ci < N - 1; ++ci)
      os << std::setw(10) << mat(ri, ci) << Char(',');
    os << std::setw(10) << mat(ri, N - 1) << Char(',') << std::endl;
  }
  return os;
}

以下のようなコードで

int main() {
  dynamic_matrix<double> mat({
      {1, 2, 3}, //
      {4, 5, 6}, //
      {7, 8, 9}  //
  });
  std::cout << mat << std::endl;
}

以下のような出力となる。

    +1.000,    +2.000,    +3.000,
    +4.000,    +5.000,    +6.000,
    +7.000,    +8.000,    +9.000,

一旦のまとめ

これ以上は長くなるので次の記事に。
クラスは以下のようになります。


#include <concepts>
#include <memory>
#include <cassert>
#include <iostream>
#include <iomanip>

template <typename T>
concept matrix_type = requires(T &a, T &b) {
  { a + b } -> std::convertible_to<T>;
  { a - b } -> std::convertible_to<T>;
  { a * b } -> std::convertible_to<T>;
  { a / b } -> std::convertible_to<T>;
  a = b;
  a = std::move(b);
};

static_assert(matrix_type<int>);
static_assert(matrix_type<float>);
static_assert(matrix_type<double>);

template <matrix_type T, typename Alloc = std::allocator<T>>
class dynamic_matrix {
public:
  using allocator_type  = std::allocator_traits<Alloc>::template rebind_alloc<T>;
  using value_type      = std::allocator_traits<allocator_type>::value_type;
  using pointer         = std::allocator_traits<allocator_type>::pointer;
  using const_pointer   = std::allocator_traits<allocator_type>::const_pointer;
  using reference       = value_type &;
  using const_reference = const value_type &;
  using size_type       = std::allocator_traits<allocator_type>::size_type;
  using index_type      = size_type;
  //dynamic_matrixのコンストラクタ
  template <std::convertible_to<T> U>
  constexpr dynamic_matrix(std::initializer_list<std::initializer_list<U>> mat,
                           const allocator_type &alloc = allocator_type())
      : //allocatorに代入演算子が存在しない事もあるので初期化はここでのみ。
      m_alloc(alloc) {
    size_type M = mat.size();
    size_type N = mat.begin()->size();
    resize(M, N);
    index_type ri = 0;
    for (auto &&row : mat) {
      index_type ci = 0;
      for (auto &&value : row)
        (*this)(ri, ci++) = static_cast<T>(value);
      ++ri;
    }
  }
  //メモリアクセス
  constexpr reference operator()(index_type i, index_type j) &noexcept {
    assert(m_elem != nullptr);
    assert(i < m_rows);
    assert(j < m_cols);
    return m_elem[m_colcap * i + j];
  }
  constexpr const_reference operator()(index_type i, index_type j) const &noexcept {
    assert(m_elem != nullptr);
    assert(i < m_rows);
    assert(j < m_cols);
    return m_elem[m_colcap * i + j];
  }
  constexpr value_type operator()(index_type i, index_type j) const &&noexcept {
    assert(m_elem != nullptr);
    assert(i < m_rows);
    assert(j < m_cols);
    return m_elem[m_colcap * i + j];
  }
  //実際に使うサイズ
  constexpr size_type rows() const noexcept {
    return m_rows;
  }
  constexpr size_type columns() const noexcept {
    return m_cols;
  }
  constexpr size_type size() const noexcept {
    return rows() * columns();
  }
  //余分に確保しているメモリサイズ
  constexpr size_type row_capacity() const noexcept {
    return m_rowcap;
  }
  constexpr size_type column_capacity() const noexcept {
    return m_colcap;
  }
  constexpr size_type capacity() const noexcept {
    return row_capacity() * column_capacity();
  }
  //アロケータを返す
  constexpr allocator_type get_allocator() const noexcept {
    return m_alloc;
  }
  //メモリ確保
  constexpr void reserve(
      size_type rows, size_type cols, size_type rowcap, size_type colcap) {
    //rowcap,colcapはそれぞれrows,colsより小さくてはいけない
    if (rows > rowcap)
      rowcap = rows;
    if (cols > colcap)
      colcap = cols;
    //現在のcapacityより大きくメモリを確保したいときのみメモリ確保を行う
    if (row_capacity() < rowcap || column_capacity() < colcap) {
      //一時的なメモリを確保
      pointer tmp = traits::allocate(m_alloc, rowcap * colcap);
      //constexpr関数内で使う為にメモリ全て初期化する。
      for (pointer sptr = tmp, eptr = tmp + rowcap * colcap; sptr < eptr; ++sptr)
        traits::construct(m_alloc, sptr);
      //すでにm_elemでメモリを確保されている場合にのみtmpにデータをコピー
      if (m_elem) {
        size_type dst_rows = std::min(rows, m_rows);
        size_type dst_cols = std::min(cols, m_cols);

        for (size_type ri = 0; ri < dst_rows; ++ri)
          for (size_type ci = 0; ci < dst_cols; ++ci)
            tmp[colcap * ri + ci] = (*this)(ri, ci);
        traits::deallocate(m_alloc, m_elem, capacity());
      }
      m_elem = tmp;
      m_rowcap = rowcap;
      m_colcap = colcap;
    }

    m_rows   = rows;
    m_cols   = cols;
  }
  //容量のみの変更を想定
  constexpr void reserve(size_type rowcap, size_type colcap) {
    reserve(m_rows, m_cols, rowcap, colcap);
  }
  //サイズの変更を想定。
  //ただし現在の容量を超えるサイズであればメモリを再確保する。
  constexpr void resize(size_type rows, size_type cols) {
    reserve(rows, cols, m_rowcap, m_colcap);
  }
  //デストラクタもconstexprとする事。
  constexpr ~dynamic_matrix() {
    deallocate();
  }
  //rows()*columns()の範囲を指定の値で初期化
  constexpr void fill(const T& value) {
    for (size_type ri = 0; ri < m_rows; ++ri)
      for (size_type ci = 0; ci < m_cols; ++ci)
        (*this)(ri, ci) = value;
  }
private:
  constexpr void deallocate(bool post_processing = true) {
    //重要: m_elemがnullptrの場合、allocateされていないか、
    //すでにdeallocateされているものとする。
    if (m_elem) {
      //後処理
      for (pointer sptr = m_elem, eptr = m_elem + capacity(); sptr < eptr; ++sptr)
        traits::destroy(m_alloc, sptr);
      //メモリ解放
      traits::deallocate(m_alloc, m_elem, capacity());
      m_elem = nullptr;
    }
    //後始末をしないほうがいいときがあるのでpost_processingで分岐する。
    //基本は後始末をしたほうがいいのでデフォルトではtrueである。
    if (post_processing) {
      m_elem   = nullptr;
      m_rows   = 0;
      m_cols   = 0;
      m_rowcap = 0;
      m_colcap = 0;
    }
  }

private:
  using traits = std::allocator_traits<allocator_type>;
  //m_elemがnullptrであるときallocatorされていない、
  //もしくはdeallocate済みで有ることを保証すること。
  //関数から抜けたときにm_elemがnullptr以外であれば
  //allocateされているという事も保証する事
  //constexpr対応の為に必須の要件とする。
  pointer m_elem         = nullptr;
  allocator_type m_alloc = {};
  size_type m_rows       = 0;
  size_type m_cols       = 0;
  size_type m_rowcap     = 0;
  size_type m_colcap     = 0;
};

template <typename Char, matrix_type T, typename Alloc>
std::basic_ostream<Char> &operator<<(
  std::basic_ostream<Char> &os, const dynamic_matrix<T, Alloc> &mat) {
  using size_type = dynamic_matrix<T, Alloc>::size_type;
  size_type M     = mat.rows();
  size_type N     = mat.columns();
  os << std::showpos << std::fixed << std::setprecision(3);
  for (size_type ri = 0; ri < M; ++ri) {
    for (size_type ci = 0; ci < N - 1; ++ci)
      os << std::setw(10) << mat(ri, ci) << Char(',');
    os << std::setw(10) << mat(ri, N - 1) << Char(',') << std::endl;
  }
  return os;
}
補足

現時点でconstexpr関数内で使用は可能である。

constexpr double test() {
  dynamic_matrix<double> mat({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}});
  return mat(1, 1);
}
int main(){
  constexpr auto res = test();//res: 5
  //以下のようなconstexpr変数はC++20時点でも不可である。
  //constexpr dynamic_matrix<double> mat1({{1,2,3},{4,5,6}});
}

最終的なまとめは以下のリポジトリ参照

https://github.com/ygsiro/portal

Discussion

ログインするとコメントできます