A single header file template container class of pseudo 2 dimensions row-major and column-major matrix

C++ における行列コンテナ sstd::mat_c<T>/sstd::mat_r<T> の提案

Abstract

　C++ において，画像処理や線形代数演算を行う場合，動的確保された固定長 2 次元配列クラスが用いられる．一般に，画像では，行方向にデータが格納されているため，row-major (行優先) の配列が用いられる．また，数値行列では，固有値計算における固有ベクトルなど，必要な計算結果が列方向に現れるため，column-major (列優先) の配列が用いられる．例えば，画像処理ライブラリである OpenCV の cv::Mat クラスは row-major として，線形代数演算ライブラリである Eigen の Matrix クラスは column-major として※1，それぞれ設計されている．また，可変長 2 次元配列を動的確保する場合，std::vector<std::vector<T>> が用いられる．これは，要素の追加や削除が発生する場合に有効である．
　このように，C++ において 2 次元配列を利用する場合は，用途に応じて型を使い分けることが望ましい．しかし，実態として，OpenCV や Eigen を用いない場合は，殆ど std::vector<std::vector<T>> が用いられている．これは，現状の STL が標準的な行列コンテナとして std::matrix<T> を提供しておらず，また，代替クラスの設計にはコストが掛かることに起因する．基本的に，std::vector<std::vector<T>> は row-major であるため，column-major のデータを扱うには不向きである．加えて，各行ベクトルは独立してメモリ確保されるため※2，メモリアドレスは不連続となる．これらは，キャッシュミスを誘発する．更に，例え固定長であっても，std::vector<T> は，現在確保しているメモリ長と挿入されている要素数を保持するため，余分にメモリを消費する．これらの問題を解消するため，本投稿では，標準的な行列コンテナとして，sstd::mat_r<T> (row-major: r) と sstd::mat_c<T> (column-major: c) を提案する．

※1. Eigen は，行列サイズが小さい場合，静的確保を用いる ^[1]．
※2. 行ベクトルのコンストラクタは，それぞれの行ごとに別に呼び出される．

Introduction

  CPU がデータを読み込むとき，読み込み対象のデータだけでなく，空間的局所性 ^[2] を期待して，隣接するデータも一緒にキャッシュメモリへ読み込まれる．これは，CPU の処理速度に対して，メモリからのデータ読み込みが非常に遅いため，より CPU に近く高速なメモリに，使用頻度の高そうなデータを置いておくほうが，効率がよいからである．このため，キャッシュミスの少ないソフトウェアを開発するには，実際に使用するデータが連続するように配置し，空間局所性を高めることが望ましい．このように，データの連続性を考えるとき，小さなメモリ領域を複数回に分けて動的確保することは，メモリ領域の連続性が保証されないため，望ましくない．即ちに，リンクリストのような，動的に確保されたメモリ領域を示すポインタは，人為的に調整しない限り連続性が保証されず，保証されない場合，ポインタを辿る度にキャッシュミスが発生する恐れがある．したがって，ソフトウェア開発において，データ領域を連続したメモリ空間上に一括して確保することは，キャッシュミスを削減しようとする試みにおいて，一般的である．これには，安易なポインタ接続を削減する必要がある．

Previous method

　C++ において 2 次元配列を確保する方法の一つは，

#include <iostream>
#include <vector>

template<typename T>
std::vector<std::vector<T>> vvalloc(const int rows, const int cols){
    return std::vector<std::vector<T>>(rows, std::vector<T>(cols));
}

template<typename T>
void vvprint(const std::vector<std::vector<T>>& rhs){
    for(int p=0; p<rhs.size(); p++){
        std::cout << "[ ";
        for(int q=0; q<rhs[p].size(); q++){
            std::cout << rhs[p][q] << ' ';
        }
        std::cout << "]" << std::endl;
    }
}

int main(){
    std::vector<std::vector<int>> vvec = vvalloc<int>(3, 3);
    
    int i=0;
    for(int p=0; p<vvec.size(); p++){
        for(int q=0; q<vvec[p].size(); q++){
            vvec[p][q] = i++;
        }
    }
    
    vvprint<int>(vvec);
    return 0;
}

実行結果

$ g++ main.cpp
$ ./a.out
[ 0 1 2 ]
[ 3 4 5 ]
[ 6 7 8 ]

Run on Ideone.com

である．このように，std::vector<T> の重ねがけによる実装は，簡便なため多用される．このとき生成されるのとほぼ同じデータ構造を C 言語で再現すると，

#include <stdio.h>
#include <stdlib.h>

struct vec_int{
    int allocate_size;
    int size; // used size
    int* val;
};
void valloc(struct vec_int* p_vec, int size){
    p_vec->allocate_size = size;
    p_vec->size = size;
    p_vec->val = (int*)malloc(size*sizeof(int));
    return p_vec;
}

struct vvec_int{
    int allocate_size;
    int size; // used size
    struct vec_int* val;
};
void vvalloc(struct vvec_int* p_vvec, int rows, int cols){
    p_vvec->allocate_size=rows;
    p_vvec->size=rows;
    p_vvec->val = (struct vec_int*)malloc(rows*sizeof(struct vec_int));
    for(int i=0; i<rows; i++){
        valloc(&p_vvec->val[i], cols);
    }
    return p_vvec;
}

void vvprint(const struct vvec_int* p_vvec){
    for(int p=0; p<p_vvec->size; p++){
        printf("[ ");
        for(int q=0; q<p_vvec->val[p].size; q++){
            printf("%d ", p_vvec->val[p].val[q]);
        }
        printf("]\n");
    }
}

int main(){
    struct vvec_int vvec;
    vvalloc(&vvec, 3, 3);
    
    int i=0;
    for(int p=0; p<vvec.size; p++){
        for(int q=0; q<vvec.val[p].size; q++){
            vvec.val[p].val[q] = i++;
        }
    }
    
    vvprint(&vvec);
    return 0;
}

実行結果

$ gcc main.cpp
$ ./a.out
[ 0 1 2 ]
[ 3 4 5 ]
[ 6 7 8 ]

Run on Ideone.com

となる (ただし，メモリの解放は省略．また，あくまで std::vector<std::vector<int>> のデータ構造の再現であり，C 言語実装としては不的確)．このようなデータ構造は，固定長 2 次元配列 (行列型) を必要とする場合において過剰である．また，メモリ空間の連続性の観点からも，好ましくない．

Purpose

  std::vector<T> の重ねがけによる多次元配列の確保は，要素の追加や削除が必要な応用において，有用である．しかしながら，行列のように，サイズが固定される用途においては，メモリ空間の不連続性が際立つ．また，確保済みメモリサイズと，使用済みメモリサイズを，行ごとに保存する分，不必要にメモリを消費する．これらの問題を解決するため，C 言語で再現するデータ構造が

#include <stdio.h>
#include <stdlib.h>

struct mat_int{
    int rows;
    int cols; // used size
    int size;
    int* val;
};
void mat_alloc(struct mat_int* p_mat, int row_size, int col_size){
    p_mat->rows = row_size;
    p_mat->cols = col_size;
    p_mat->size = row_size * col_size;
    p_mat->val = (int*)malloc(p_mat->size * sizeof(int));
    return p_mat;
}

void mat_print(const struct mat_int* p_mat){
    for(int p=0; p<p_mat->rows; p++){
        printf("[ ");
        for(int q=0; q<p_mat->cols; q++){
            printf("%d ", p_mat->val[p_mat->cols*p + q]);
        }
        printf("]\n");
    }
}

int main(){
    struct mat_int mat;
    mat_alloc(&mat, 3, 3);
    
    int i=0;
    for(int i=0; i<mat.size; i++){
        mat.val[i] = i;
    }
    
    mat_print(&mat);
    return 0;
}

実行結果

$ gcc main.cpp
$ ./a.out
[ 0 1 2 ]
[ 3 4 5 ]
[ 6 7 8 ]

Run on Ideone.com

となる C++ のポータブルな行列コンテナを実現したい．このとき，コンテナのインターフェースは，

#include <iostream>
#include "sstd_mat_r.hpp"

template<typename T>
void mat_print(const sstd::mat_r<T>& rhs){
    for(int p=0; p<rhs.rows(); p++){
        std::cout << "[ ";
        for(int q=0; q<rhs.cols(); q++){
            std::cout << rhs(p, q) << ' ';
        }
        std::cout << "]" << std::endl;
    }
}

int main(){
    sstd::mat_r<int> mat = {{0, 1, 2}, {3, 4, 5}, {6, 7, 8}};

    mat_print(mat);

    return 0;
}

となるように設計する．

Required specifications

　本投稿における要求仕様の概略は，
・() 演算子をオーバーロードして，1 次元配列へのアクセスを，擬似的に 2 次元とする．
・{} を用いたコンテナ型の初期化．
・自然なコンテナの初期化とコピー．
・名前空間による名称衝突の回避．
である．
　特に，コンテナの初期化やコピーは，複数通りあるため，注意が必要である．

Implementation

  要求に沿った実装を下記に示す．本実装は，拙ライブラ sstd より過剰な実装を削除し，シングルヘッダライブラリとした．また，実装はコメントを含めて 155 行である．なお，column-major 版の実装である sstd_mat_c.hpp の掲載は省略するため，必要であれば後述の GitHub のリンクを辿って頂きたい．

./sstd_mat_r.hpp

#pragma once
#include <initializer_list>

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

typedef unsigned int uint;

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

namespace sstd{
    template <typename T> class mat_r;
    template <typename T> void copy(sstd::mat_r<T>& lhs, const sstd::mat_r<T>& rhs);
    template <typename T> void move(sstd::mat_r<T>& lhs,       sstd::mat_r<T>& rhs);
    template <typename T> void swap(sstd::mat_r<T>& lhs,       sstd::mat_r<T>& rhs);
    
    template <typename T> sstd::mat_r<T> Tr       (const sstd::mat_r<T>& rhs); // lhs = Transpose(rhs)
    template <typename T> void           Tr_myself(      sstd::mat_r<T>& rhs); // Transpose(rhs)
}

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

// type martix (row major)
template <typename T>
class sstd::mat_r{
private:
    T* _pMatT;
    uint _rows; // row size
    uint _cols; // column size
    uint _size; // number of elements

public:
    inline mat_r(){ _rows=0; _cols=0; _size=0; _pMatT=0; }
    inline mat_r(const std::initializer_list<std::initializer_list<T>>& rhs){ // called by "sstd::mat_r<T> mat = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};".
        _rows = rhs.size();
        if(_rows==0){ _cols=0; _size=0; return; }
        _cols  = (rhs.begin()[0]).size();
        _size  = _rows*_cols;
        _pMatT = new T[_size];
        
        const std::initializer_list<T>* pRhs=rhs.begin();
        for(uint p=0; p<_rows; p++){
            const T* ppRhs=pRhs[p].begin();
            for(uint q=0; q<_cols; q++){
                if(q>=pRhs[p].size()){break;}
                _pMatT[p*_cols + q]=ppRhs[q];
            }
        }
    }
    inline mat_r(const class mat_r&  rhs){ _rows=0; _cols=0; _size=0; _pMatT=0; sstd::copy<T>(*this, rhs); } // called by "sstd::mat_r<T> buf1(N, N); sstd::mat_r<T> buf2(buf1);"
    inline mat_r(      class mat_r&& rhs){ _rows=0; _cols=0; _size=0; _pMatT=0; sstd::move<T>(*this, rhs); } // called by "return std::move(rhs);" or "std::swap(buf1, buf2)".
    inline mat_r(const uint row_size, const uint col_size){
        _rows = row_size;
        _cols = col_size;
        _size = _rows * _cols;
        _pMatT = new T[_size];
        if(_pMatT==0){ printf("ERROR: sstd::mat_r<T>: This pointer is not allocated.\n"); exit(-1); }
    }
    inline ~mat_r(){ delete[] _pMatT; }
    
    // R: read only
    inline const uint rows() const { return _rows; }
    inline const uint cols() const { return _cols; }
    inline const uint size() const { return _size; }
    
    // RW: read and write
    inline T*& pMatT_RW(){ return _pMatT; }
    inline uint& rows_RW(){ return _rows; }
    inline uint& cols_RW(){ return _cols; }
    inline uint& size_RW(){ return _size; }
    
    class mat_r& operator=(const class mat_r& rhs){ sstd::copy<T>(*this, rhs); return *this; } // called by "lhs = sstd::mat_r<T>(3, 3);".
    
    inline       T& operator[](const uint i)       { return _pMatT[      i]; }
    inline       T& operator[](const  int i)       { return _pMatT[(uint)i]; }
    inline const T& operator[](const uint i) const { return _pMatT[      i]; }
    inline const T& operator[](const  int i) const { return _pMatT[(uint)i]; }

    inline       T& operator()(const uint p, const uint q)       { return _pMatT[_cols*      p +       q]; }
    inline       T& operator()(const  int p, const uint q)       { return _pMatT[_cols*(uint)p +       q]; }
    inline       T& operator()(const uint p, const  int q)       { return _pMatT[_cols*      p + (uint)q]; }
    inline       T& operator()(const  int p, const  int q)       { return _pMatT[_cols*(uint)p + (uint)q]; }
    inline const T& operator()(const uint p, const uint q) const { return _pMatT[_cols*      p +       q]; }
    inline const T& operator()(const  int p, const uint q) const { return _pMatT[_cols*(uint)p +       q]; }
    inline const T& operator()(const uint p, const  int q) const { return _pMatT[_cols*      p + (uint)q]; }
    inline const T& operator()(const  int p, const  int q) const { return _pMatT[_cols*(uint)p + (uint)q]; }
};

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

template <typename T>
inline void copy_withoutAllocation(sstd::mat_r<T>& lhs, const sstd::mat_r<T>& rhs){
    for(uint p=0; p<rhs.rows(); p++){
        for(uint q=0; q<rhs.cols(); q++){
            lhs(p, q) = rhs(p, q);
        }
    }
}
template <typename T>
inline void sstd::copy(sstd::mat_r<T>& lhs, const sstd::mat_r<T>& rhs){

    if(lhs.size() != rhs.size()){
        delete[] lhs.pMatT_RW(); lhs.pMatT_RW() = 0;
        if(rhs.size()!=0){ lhs.pMatT_RW() = new T[rhs.size()]; }
    }

    lhs.rows_RW() = rhs.rows();
    lhs.cols_RW() = rhs.cols();
    lhs.size_RW() = rhs.size();

    copy_withoutAllocation<T>(lhs, rhs);
}

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

template <typename T>
inline void sstd::move(sstd::mat_r<T>& lhs, sstd::mat_r<T>& rhs){    
    lhs.rows_RW() = rhs.rows(); rhs.rows_RW() = 0;
    lhs.cols_RW() = rhs.cols(); rhs.cols_RW() = 0;
    lhs.size_RW() = rhs.size(); rhs.size_RW() = 0;

    delete[] lhs.pMatT_RW();
    lhs.pMatT_RW() = rhs.pMatT_RW();
    rhs.pMatT_RW() = 0;
}

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

template <typename T>
void swap(sstd::mat_r<T>& lhs, sstd::mat_r<T>& rhs){
    uint rowsBuf=lhs.rows(); lhs.rows_RW()=rhs.rows(); rhs.rows_RW()=rowsBuf;
    uint colsBuf=lhs.cols(); lhs.cols_RW()=rhs.cols(); rhs.cols_RW()=colsBuf;
    uint sizeBuf=lhs.size(); lhs.size_RW()=rhs.size(); rhs.size_RW()=sizeBuf;
    
    T* pMatBuf=lhs.pMat(); lhs.pMatT_RW()=rhs.pMat(); rhs.pMatT_RW()=pMatBuf;
}

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

template <typename T>
sstd::mat_r<T> sstd::Tr(const sstd::mat_r<T>& rhs){
    sstd::mat_r<T> lhs(rhs.cols(), rhs.rows());
    for(uint p=0; p<rhs.rows(); p++){
        for(uint q=0; q<rhs.cols(); q++){
            lhs(q, p) = rhs(p, q);
        }
    }
    return lhs;
}
template <typename T>
void sstd::Tr_myself(sstd::mat_r<T>& rhs){
    sstd::mat_r<T> lhs = sstd::Tr(rhs);
    sstd::move(rhs, lhs);
}

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

　上記の実装について，順に説明する．
　まず，コンストラクタとディストラクタについて説明する．sstd::mat_r.hpp のコンストラクタは，

    inline mat_r(){ _rows=0; _cols=0; _size=0; _pMatT=0; }
    inline mat_r(const std::initializer_list<std::initializer_list<T>>& rhs){ // called by "sstd::mat_r<T> mat = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};".
        _rows = rhs.size();
        if(_rows==0){ _cols=0; _size=0; return; }
        _cols  = (rhs.begin()[0]).size();
        _size  = _rows*_cols;
        _pMatT = new T[_size];
        
        const std::initializer_list<T>* pRhs=rhs.begin();
        for(uint p=0; p<_rows; p++){
            const T* ppRhs=pRhs[p].begin();
            for(uint q=0; q<_cols; q++){
                if(q>=pRhs[p].size()){break;}
                _pMatT[p*_cols + q]=ppRhs[q];
            }
        }
    }
    inline mat_r(const class mat_r&  rhs){ _rows=0; _cols=0; _size=0; _pMatT=0; sstd::copy<T>(*this, rhs); } // called by "sstd::mat_r<T> buf1(N, N); sstd::mat_r<T> buf2(buf1);"
    inline mat_r(      class mat_r&& rhs){ _rows=0; _cols=0; _size=0; _pMatT=0; sstd::move<T>(*this, rhs); } // called by "return std::move(rhs);" or "std::swap(buf1, buf2)".
    inline mat_r(const uint row_size, const uint col_size){
        _rows = row_size;
        _cols = col_size;
        _size = _rows * _cols;
        _pMatT = new T[_size];
        if(_pMatT==0){ printf("ERROR: sstd::mat_r<T>: This pointer is not allocated.\n"); exit(-1); }
    }

の 5 種類が実装されている．ここでは，上から順に，コンストラクタ 1., 2., ... 5. と仮称し，呼び出しタイミングを下記に示す．

コンストラクタ 1.

    inline mat_r(){ _rows=0; _cols=0; _size=0; _pMatT=0; }

の呼び出しタイミングは，

sstd::mat_r<T> mat; // ただし，変数を使用しない場合は最適化により実体化されない．
sstd::mat_r<T> mat = sstd::mat_r<T>(); // ただし，変数を使用しない場合は最適化により実体化されない．

の 2 通りである．このとき，= 演算子は呼び出されない．

コンストラクタ 2.

    inline mat_r(const std::initializer_list<std::initializer_list<T>>& rhs){ // called by "sstd::mat_r<T> mat = {{1, 2, 3}, {4, 5, 6}, {7, 8, 

の呼び出しタイミングは，

sstd::mat_r<T> mat={{0,1,2}, {3,4,5}, {6,7,8}}; // ただし，この場合 T は数値型

である．このとき，= 演算子は呼び出されない．

コンストラクタ 3.

    inline mat_r(const class mat_r&  rhs){ _rows=0; _cols=0; _size=0; _pMatT=0; sstd::copy<T>(*this, rhs); } // called by "sstd::mat_r<T> buf1(N, N);

の呼び出しタイミングは，

sstd::mat_r<T> mat(mat_dummy);
sstd::mat_r<T> mat = mat_dummy;
sstd::mat_r<T> mat = sstd::mat_r<T>(mat_dummy);

の 3 通りである．いずれも，= 演算子は呼び出されない．また，mat_dummy は，初期化済みの sstd::mat_r<T> 型変数を表す．

コンストラクタ 4.

    inline mat_r(      class mat_r&& rhs){ _rows=0; _cols=0; _size=0; _pMatT=0; sstd::move<T>(*this, rhs); } // called by "return std::move(rhs);" or "std::swap(buf1, buf2)".

の呼び出しタイミングは，

sstd::mat_r<T> mat(std::move(mat_dummy));
sstd::mat_r<T> mat = std::move(mat_dummy);
sstd::mat_r<T> mat = sstd::mat_r<T>(std::move(mat_dummy));

の 3 通りである．いずれも，= 演算子は呼び出されない．また，mat_dummy は，初期化済みの sstd::mat_r<T> 型変数を表す．

コンストラクタ 5.

    inline mat_r(const uint row_size, const uint col_size){

の呼び出しタイミングは，

sstd::mat_r<T> mat(3, 3);

のである．

　なお，関数の戻り値として

sstd::mat_r<T> mat = function_of_return_mat();

のように初期化される場合は，関数内で発生するコンストラクタ以外は，何も呼び出されない．コンストラクタ 4. さえ呼び出されないことから，コンパイラが完全に最適化していることがわかる．なお，この結果は，gcc 5.4.0 上で，単純なサンプルプログラムに対してコンパイルした結果である．

　また，ディストラクタは，

    inline ~mat_r(){ delete[] _pMatT; }

となっている．ディストラクタの呼び出しタイミングは，C++ の仕様に従う．オブジェクトを静的に扱う場合，オブジェクトの寿命は，{ から } までであるから，} を抜けたタイミングで，ディストラクタが呼び出される．これは，for 文や while 文であっても，{ から } までである．このようなメモリ管理は RAII と呼ばれる．また，動的に確保した場合，すなわち，new を用いて確保した場合は，delete を呼び出したタイミングでディストラクタが呼び出される．

　注意点として，配列は new 確保する必要がある．当然ではあるが，T にオブジェクト型を用いる場合，malloc によるメモリ確保は，コンストラクタを呼び出さない．

　次に，内部変数である配列サイズを，外部から取得するためのメンバ関数について，説明する．
　内部変数を取得するメンバ関数は，それぞれ，

    // R: read only
    inline const uint rows() const { return _rows; }
    inline const uint cols() const { return _cols; }
    inline const uint size() const { return _size; }

のように実装されている．命名は，Eigen に倣い，行数を rows()，列数を cols() と定義している．また，要素数は STL に倣い size() とした．これらの値は，基本的に，外部から変更することはないため，const メンバ関数としている．C++ では「const オブジェクトは 非 const メンバ関数を呼び出 ^[3]」せないため，注意が必要である．なお，戻り値の const は明示的に記述しているだけであり，取り除いても動作に変化はない．
　さて，"基本的には" これらの内部変数を外部から変更することはありえないが，実装を簡潔に保つため，

    // RW: read and write
    inline T*& pMatT_RW(){ return _pMatT; }
    inline uint& rows_RW(){ return _rows; }
    inline uint& cols_RW(){ return _cols; }
    inline uint& size_RW(){ return _size; }

のように，外部からアクセスするためのメンバ関数を用意する．ここでは，値を変更する必要があるため，戻り値には参照渡しを用いる．

　オブジェクトのコピーを行う場合には，= 演算子も定義されている必要がある．ここでは，= 演算子は，

    class mat_r& operator=(const class mat_r& rhs){ sstd::copy<T>(*this, rhs); return *this; } // called by "lhs = sstd::mat_r<T>(3, 3);".

のように定義されている．呼び出しタイミングは，

sstd::mat_r<T> matR={{0,1,2}, {3,4,5}, {6,7,8}}; // コンストラクタ 2. の呼び出し．(ただし，この場合 T は数値型)
sstd::mat_r<T> mat; // コンストラクタ 1. の呼び出し．
mat = matR; // = 演算子の呼び出し．

である．

　内部配列への通常アクセスには，[] 演算子を，

    inline       T& operator[](const uint i)       { return _pMatT[      i]; }
    inline       T& operator[](const  int i)       { return _pMatT[(uint)i]; }
    inline const T& operator[](const uint i) const { return _pMatT[      i]; }
    inline const T& operator[](const  int i) const { return _pMatT[(uint)i]; }

のように定義して対処する．sstd::mat_r<T> の T 型は，built-in 型だけでなく，任意のオブジェクト型に対して定義されるため，戻り値は参照渡しである必要がある．これは，T 型を巨大なイブジェクト型としてインスタンス化した場合に，多量のメモリコピーが発生することを防ぐためである．参照渡しは，コンパイル時にポインタとして処理されるため，このコストは，ポインタのコピーと同程度と考えることができる (実装依存のため，コンパイラ依存とも考えられる)．ポインタの大きさはアドレスサイズと等価であるから，64 bits CPU では 64 bits である．また，const オブジェクトとして用いられる場合に備えて，const メンバ関数を用意する必要がある．

　擬似 2 次元配列としてのアクセスには，() 演算子を，

    inline       T& operator()(const uint p, const uint q)       { return _pMatT[_cols*      p +       q]; }
    inline       T& operator()(const  int p, const uint q)       { return _pMatT[_cols*(uint)p +       q]; }
    inline       T& operator()(const uint p, const  int q)       { return _pMatT[_cols*      p + (uint)q]; }
    inline       T& operator()(const  int p, const  int q)       { return _pMatT[_cols*(uint)p + (uint)q]; }
    inline const T& operator()(const uint p, const uint q) const { return _pMatT[_cols*      p +       q]; }
    inline const T& operator()(const  int p, const uint q) const { return _pMatT[_cols*(uint)p +       q]; }
    inline const T& operator()(const uint p, const  int q) const { return _pMatT[_cols*      p + (uint)q]; }
    inline const T& operator()(const  int p, const  int q) const { return _pMatT[_cols*(uint)p + (uint)q]; }

のように定義して対処する．上記と同様に，配列の書き換え用 non-const メンバ関数と，const オブジェクトとして用いられる場合に備えて const メンバ関数を，それぞれ用意する必要がある．

Results

　サンプルコードと，その実行結果を示す．

main.cpp

#include <iostream>
#include "./sstd_mat_c.hpp"
#include "./sstd_mat_r.hpp"

typedef unsigned int uint;

template<typename T>
void mat_print(const sstd::mat_c<T>& rhs){
    for(uint p=0; p<rhs.rows(); p++){
        std::cout << "[ ";
        for(uint q=0; q<rhs.cols(); q++){
            std::cout << rhs(p, q) << ' ';
        }
        std::cout << "]" << std::endl;
    }
    std::cout << std::endl;
}

template<typename T>
void mat_print(const sstd::mat_r<T>& rhs){
    for(uint p=0; p<rhs.rows(); p++){
        std::cout << "[ ";
        for(uint q=0; q<rhs.cols(); q++){
            std::cout << rhs(p, q) << ' ';
        }
        std::cout << "]" << std::endl;
    }
    std::cout << std::endl;
}

int main(){
    // TEST of sstd::mat_c<int>.
    {
        sstd::mat_c<int> matC2 = {{0, 1, 2}, {3, 4, 5}, {6, 7, 8}};
        sstd::mat_c<int> matC1 = matC2;
        sstd::mat_c<int> matC3(matC1);
        sstd::mat_c<int> matC4(std::move(matC3));
        sstd::mat_c<int> matC5(3, 3);
        for(uint i=0; i<matC5.size(); i++){ matC5[i]=matC4[i]; }
        sstd::mat_c<int> matC5_2(3, 3);
        for(uint q=0; q<matC5_2.cols(); q++){
            for(uint p=0; p<matC5_2.rows(); p++){
                matC5_2(p, q) = matC5(p, q);
            }
        }
        mat_print(matC5_2);
    }
    
    // TEST of sstd::mat_r<int>.
    {
        sstd::mat_r<int> matR2 = {{0, 1, 2}, {3, 4, 5}, {6, 7, 8}};
        sstd::mat_r<int> matR1 = matR2;
        sstd::mat_r<int> matR3(matR1);
        sstd::mat_r<int> matR4(std::move(matR3));
        sstd::mat_r<int> matR5(3, 3);
        for(uint i=0; i<matR5.size(); i++){ matR5[i]=matR4[i]; }
        sstd::mat_r<int> matR5_2(3, 3);
        for(uint p=0; p<matR5_2.rows(); p++){
            for(uint q=0; q<matR5_2.cols(); q++){
                matR5_2(p, q) = matR5(p, q);
            }
        }
        mat_print(matR5_2);
    }
    return 0;
}

実行結果

$ make
$ ./a.out 
[ 0 1 2 ]
[ 3 4 5 ]
[ 6 7 8 ]

[ 0 1 2 ]
[ 3 4 5 ]
[ 6 7 8 ]

Conclusion

  拙ライブラリ sstd より，余剰機能を削除した，任意型 T に対する row-major および column-major 擬似 2 次元配列コンテナのポータブル版を実装した．本ヘッダライブラリは，オリジナル版から行列に対する算術演算子の定義や，要素演算の定義が削除し，単一ヘッダライブラリとして，最小限の機能を実装している．

Source code

本投稿で提示したコードの URL を下記に示す．
github.com/admiswalker/blog_2019_04_matCR

拙ライブラリ sstd の URL を示す．今後変更される場合はこちらが変更される．
github.com/admiswalker/SubStandardLibrary-SSTD-

References

  [1] C++行列計算ライブラリEigen入門 - Quitta - 2019年04月05日閲覧
  [2] キャッシュ (コンピュータシステム) - Wikipedia - 2019年04月05日閲覧
  [3] constメンバ関数 - WisdomSoft (旧) > プログラミング > C++入門 - 2019年04月08日閲覧