数値計算で慣れようC++　Part2. 行列積を実装しよう #C++

この記事は CAMPHOR- Advent Calendar 2023 16日目の記事です。

前回までのあらすじ

C++で有理数型を実装した!
今回は行列積を計算するコードを書いてC++のメモリ管理やテンプレートについて知ろう!

※前回の記事を読まなくてもなんとかなります。

行列積を実装しよう

行列積の実装を通じて、さらにC++の機能をみていきましょう。
なお、行/列ベクトルは単なる配列として定義し、行列は行ベクトルのポインタの配列として定義します。

C言語の実装

まずはC言語で実装してみます。
もちろんこのコードはC++用のコンパイラでもコンパイルできます。

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

double** matrix_mul(double** m1, double** m2, int l, int m, int n){
    double** mul = (double**)calloc(l, sizeof(double*));
    for(int i = 0; i < l; i++){
        mul[i] = (double*)calloc(n, sizeof(double));
        for(int j = 0; j < n; j++){
            for(int k = 0; k < m; k++){
                mul[i][j] += m1[i][k] * m2[k][j];
            }
        }
    }
    return mul;
}

int main(){
    double A0[3] = {1,2,3};
    double A1[3] = {4,5,6};
    double A2[3] = {7,8,9};
    double* A[3] = {A0, A1, A2};
    for(int i = 0; i < 3; i++)
        printf("%lf, %lf, %lf\n", A[i][0], A[i][1], A[i][2]);
// >1.000000, 2.000000, 3.000000
// >4.000000, 5.000000, 6.000000
// >7.000000, 8.000000, 9.000000
        
    double** x = matrix_mul(A,A,3,3,3);
    for(int i = 0; i < 3; i++)
        printf("%lf, %lf, %lf\n", x[i][0], x[i][1], x[i][2]);
// >30.000000, 36.000000, 42.000000
// >66.000000, 81.000000, 96.000000
// >102.000000, 126.000000, 150.000000
        
    return 0;
    
}

m1が左の行列(l×m)で、m2が右の行列(m×n)です。
また、行列のほかに行列のサイズを引数にとっています。
内部で行っている処理は単純に3重のforループで行列積を計算しているだけです。
main関数内では適当な正方行列の二乗を計算しています。
計算結果はnumpyで適当に検算しました。

C++を用いてより良く実装する方針を２つ示します。
１つ目は、メモリの動的確保に関してです。
C言語ではmallocやcallocでメモリを確保し、freeで解放していました。
しかし、C++では代わりにnewとdeleteを用います。

2つ目はコードの可搬性です。
ここで定義されている行列積の関数はdoubleで表現された数値のみを対象にしています。
ですが、せっかくなら前回実装した有理数型に対しても同様の行列演算をしたいですよね？
C++で様々な型のデータに共通の処理を行いたい、そんなモチベーションに応えるのがテンプレートという機能です。

※前回の記事を読んでない方はfloatやbinary128(4倍精度)などに置き換えてください

newとdelete

C++では、メモリを動的に用いる際にnewやdeleteを用います。
前回定義したrational型を例に試してみましょう¹。

auto a = new rational(1,3);
// a := 1/3
std::cout << a->to_double() << std::endl;
// > 0.33333
auto b = new rational(2);
// b := 2

delete a;

auto a_ = a->to_double();
// Undefined

ご覧の通り、newというキーワードでインスタンスのためのメモリ領域を確保し、コンストラクタを呼び出します。
mallocと同様に、ポインタが帰ってきます。構造体のポインタと同様に扱えます。
また、deleteというキーワードでデストラクタを呼び出し、インスタンスに割り当てられていたメモリ領域を解放します。

以上の話を踏まえて、冒頭の行列積のコードをC++で実装すると、以下のようになります。

#include <iostream>


double** matrix_mul(double** m1, double** m2, int l, int m, int n){
    auto mul = new double*[l];
    for(int i = 0; i < l; i++){
        mul[i] = new double[n];
        for(int j = 0; j < n; j++){
            mul[i][j] = double(0);
            for(int k = 0; k < m; k++){
                mul[i][j] = mul[i][j] + m1[i][k] * m2[k][j];
            }
        }
    }
    return mul;
}

int main(){
    double A0[3] = {1,2,3};
    double A1[3] = {4,5,6};
    double A2[3] = {7,8,9};
    double* A[3] = {A0, A1, A2};
    for(int i = 0; i < 3; i++)
        std::cout << A[i][0] << " " << A[i][1] << " " << A[i][2] << std::endl;
    double** x = matrix_mul(A,A,3,3,3);
    for(int i = 0; i < 3; i++)
        std::cout << x[i][0] << " " << x[i][1] << " " << x[i][2] << std::endl;
    return 0;
    
}

newはクラスだけでなく、intやポインタなどの組み込み型や配列を定義する際にも使えます(当然ではありますが)。
また、気をつけてほしいのはクラスの配列を定義するときです。
上のコードでは、mul[i] = new double[n];とdoubleの配列を定義したあとに
mul[i][j] = double(0);と各々の要素を定義しています。
~~newでクラスの配列を定義した際は、各々の要素についてコンストラクタは呼ばれません。~~12/17追記: コード自体は動きますが記載に誤りがあります。近日中に修正致します。

ここではdouble型を用いているため各々の要素を定義する必要は必ずしもありませんが、
後にrational型を使うときにこれを忘れていると痛い目にあいます。

他にも、C++にはメモリ管理を便利にするためのスマートポインタと呼ばれる仕組みがあります。
若干ややこしいため、以下で軽く説明します。

補足　スマートポインタ

行列をポインタ配列で管理すると、範囲外アクセスやallocし忘れ、freeし忘れなどのエラーが頻発します[^2]。そういった事故を防ぐため、C++ではスマートポインタという仕組みが用意されています。以下のコードを見てみましょう。これは、スマートポインタを用いて先程の行列積のコードを書き直したものです。

#include <memory>

std::unique_ptr<std::unique_ptr<double[]>[]> matrix_mul(double** m1, double** m2, int l, int m, int n){
    auto mul = std::make_unique<std::unique_ptr<double[]>[]>(l);
    for(int i = 0; i < l; i++){
        mul[i] = std::make_unique<double[]>(n);
        for(int j = 0; j < n; j++){
            for(int k = 0; k < m; k++){
                mul[i][j] = mul[i][j] + m1[i][k] * m2[k][j];
            }
        }
    }
    return mul;
}
/*--------------------------------------------
 --------------------------------------------*/

auto x = matrix_mul(m1,m2,l,m,n);

ポインタの変わりにstd::unique_ptrという型を用いています。
ポインタ配列(ポインタのポインタ)を扱っているため、スマートポインタが入れ子状になっています。
このスマートポインタはuniqueポインタというもので、所有権を共有しない対象に用いられます。
定義の際に型や要素数を指定すると、勝手にコンストラクタを呼んで初期化してくれる優れものです。
さらに、使わなくなった変数のメモリを適切に解放してくれるため、freeし忘れも防止できます²。
さらに、変数の所有権を共有したい場合にはsharedポインタやweakポインタなども用意されています³。

テンプレート

異なる型に対して同じ処理をしたい場合、基本的にはもう一つ別に関数を定義する必要があります⁴。
でも、そんな冗長なことを書きたくないですよね。
C言語のエキスパートなら、マクロを使うかもしれませんね。
しかし、C++にはもっとクレバーな方法が用意されています。
それが、テンプレートです。
以下のコードを見てください。

#include <iostream>

template <class T>
T** matrix_mul(T** m1, T** m2, int l, int m, int n){
    auto mul = new T*[l];
    for(int i = 0; i < l; i++){
        mul[i] = new T[n];
        for(int j = 0; j < n; j++){
            mul[i][j] = T(0);
            for(int k = 0; k < m; k++){
                mul[i][j] = mul[i][j] + m1[i][k] * m2[k][j];
            }
        }
    }
    return mul;
}

int main(){
    double A0[3] = {1,2,3};
    double A1[3] = {4,5,6};
    double A2[3] = {7,8,9};
    double* A[3] = {A0, A1, A2};
    for(int i = 0; i < 3; i++)
        std::cout << A[i][0] << " " << A[i][1] << " " << A[i][2] << std::endl;
    double** x = matrix_mul<double>(A,A,3,3,3);
    for(int i = 0; i < 3; i++)
        std::cout << x[i][0] << " " << x[i][1] << " " << x[i][2] << std::endl;
    return 0;
    
}

コードの骨子は変わりませんが、templateやTなど見慣れないものが登場していますね。
C++のコードらしくなってきました。
このtemplateの記述は直後の関数やクラスを修飾するものです。
修飾される関数やクラス内でTを用いて一般化した記述を行うと、あとからTに好きな型を指定できるようになります。

実際に、main関数内でmatrix_mul<double>(略)という形でmatrix_mul関数が呼び出されています。
これによって、doubleを入出力にとる行列積が計算されています。

また、matrix_mul関数内のjに関するforの次の行をみてください。
mul[i][j] = T(0);　とあります。
この関数では、返す行列はT**型の変数として定義され、各要素もT型の変数としてアロケートされます。
この際、doubleやfloatならcallocされた状態だと0という値を保持していることになるため、とくに初期化を行わなくとも問題ありませんでした。
しかし、テンプレートを用いるときはrational型のようにコンストラクタを呼ぶ必要がある型が利用されることもあるため、このようにコンストラクタを呼ぶ実装にしてあります。

また、テンプレートでは型以外のものを受け付けることもできます。
整数を受け付ける例は散見されます。以下のコードを見てください。

#include <iostream>

template <class T, int l, int m, int n>
T** matrix_mul(T** m1, T** m2){
    auto mul = (T**)calloc(l, sizeof(T*));
    for(int i = 0; i < l; i++){
        mul[i] = (T*)calloc(n, sizeof(T));
        for(int j = 0; j < n; j++){
            mul[i][j] = T(0);
            for(int k = 0; k < m; k++){
                mul[i][j] = mul[i][j] + m1[i][k] * m2[k][j];
            }
        }
    }
    return mul;
}


int main(){
    double A0[3] = {1,2,3};
    double A1[3] = {4,5,6};
    double A2[3] = {7,8,9};
    double* A[3] = {A0, A1, A2};
    for(int i = 0; i < 3; i++)
        std::cout << A[i][0] << " " << A[i][1] << " " << A[i][2] << std::endl;
        
    double** x = matrix_mul<double, 3, 3, 3>(A,A);
    
    for(int i = 0; i < 3; i++)
        std::cout << x[i][0] << " " << x[i][1] << " " << x[i][2] << std::endl;
    return 0;   
}

行列のサイズl,m,nを引数ではなくテンプレートで受け取る形で実装しました。
この例だとありがたみがわかりにくいかもしれませんが、クラスを定義するときなどは有用だったりします(たぶん)。

ここまでみてきたテンプレートを用いたmatrix_mul関数ですが、前回のrational型に二項演算子の + , - , * , / をオーバーロードできていればrational型でも動作します⁵。
代入のための単項演算子 = も実装できているとなおよいでしょう。
ストリームに流すためのシフト演算子も実装できているとなお便利ですが、若干難点があるためいずれ解説します。

さて、次回は行列をクラスとして定義し、いろいろとメンバ関数を実装していきましょう⁶。

別にrationalである必要性はなく、前回を読んでない人は適当なクラスで試してください。 ↩
デストラクタが呼ばれる(スコープを抜けたときを含む)際に勝手にdeleteまでしてくるんですが、割とややこしいため詳細は省きます(今回の例では関数内で定義したuniqueポインタを返していますが、スコープを抜けているのになぜ返せるの？みたいな話になってくる)。 ↩
今回のテーマではあまり所有権を共有するようなデータ構造は扱わないため割愛します。良いテーマが見つかったら説明するかも。 ↩
この辺り(ジェネリクス)は言語によってかなり異なってきます。pythonのようにとりあえず動かしてみてダメならエラーを吐くような言語もあれば、Rustのようにtraitという形で管理する言語もあります。 ↩
実は2020年に策定されたC++20にはconceptという仕組みがあって、Templateで用いる型がどんな関数を実装している必要があるかを明示的に示すことができるようになっています。が、若干煩雑(演算子オーバーロードがどうなってるか筆者が理解していない)ため割愛します。 ↩
初めからクラスとして定義しておけよ、という指摘はごもっともです。が、そこそこの規模のクラスを初見のギミックを使って定義するのも大変かなと思った次第です。 ↩

数値計算で慣れようC++ Part2. 行列積を実装しよう

前回までのあらすじ

行列積を実装しよう

C言語の実装

newとdelete

テンプレート

数値計算で慣れようC++　Part2. 行列積を実装しよう