46 行列 (二重添字を持つ数列) のデータ構造について再考する

暮に講習会に出たのだけれど、そこで久しぶりに思い出すことになって、しばし再考して楽しんだ(？)話。

結局は、過去に書いたこと (「行列はどうする」特に「二つの方法の優劣」, 「行列・2次元格子上の数値データを C のプログラムでどう扱うべきか」) を思い出しただけのような気もするが。

微分方程式のシミュレーションをしていると、行列 $A=(a_{ij})$ や2次元格子上で定義された関数値 $\left\{u(x_i, y_j)\right\}$ を記録するために二重添字を持つ数列を扱う必要が生じる。 Fortran でプログラムを書く場合は、二次元配列で扱って特に問題は生じないことが多いのだが、 C 言語で書く場合はいくつか問題が生じる。

整合配列を使いたいが…

Fortran には、古くから整合配列 (conformant array) の機能があるため、副プログラムとの間のデータのやりとりに特に問題は生じないが、 C言語の場合は、そうではない。

もともと C 言語には整合配列の機能はなかった。
配列の大きさの指定に変数が使える「可変長配列 (variable length array)」という機能が、 GCC などで試験的に導入され、C99 の規格に取り込まれた (C99 については、「プログラミング言語 C の新機能」が参考になります)。これを関数引数に使うことで実質上整合配列が利用可能になる。
ところが、可変長配列の機能は C11 ではオプションにされてしまった (脱線かもしれないが、C++ では規格に入っていない)。

GCC や LLVM では、可変長配列の機能は残っているので、それを使えば良いのかもしれないが… 個人的にはかなり気持ちが悪い (規格は尊重すべき)。

私のコンピューターの師匠は、ポインターへのポインターを使いなさい、少し余分なメモリーが必要になるけれど問題にはならない、という意見だった。

そういうわけで、個人的には次のようなコードを使うことにしていた。

#include <stdlib.h> typedef double *vector; typedef vector *matrix; matrix newmatrix(int m, int n) { int i; vector ap; matrix a; if ((a = malloc(m * sizeof(vector))) == NULL) return NULL; if ((ap = malloc(m * n * sizeof(double))) == NULL) { free(a); return NULL; } for (i = 0; i < m; i++) a[i] = ap + (i * n); return a; }

(malloc() はメモリーを動的に確保する関数である。 void *malloc(size_t) というプロトタイプ宣言が、 stdlib.h に入っている。)

このやり方は確かに少し余分のメモリーが必要になるが、評価できるところは多い。

malloc() を使うので、スタック・サイズなどを気にすることなく、大きなメモリーが取れる。
$a_{ij}$ にアクセスするには、単に a[][] で良く、 2次元配列用のコードと同じで良い。

関数に引数として渡す場合も簡単である。

/* 素朴なやり方で行列の掛け算をする */
void mulmatrix(matrix c, matrix a, matrix b, int l, int m, int n)
{
  int i, j, k;
  double s;
  for (i = 0; i < l; i++)
    for (j = 0; j < n; j++) {
       s = 0.0;
       for (k = 0; k < m; k++)
         s += a[i][k] * b[k][j];
       c[i][j] = s;
    }
}

計算速度の面でも特に遅いということはない。行列の行の交換をする場合などは、ポインターの値の交換をすれば済むので有利になることもある。
古い規格の C 言語処理系でも(小さな修正で)使える。

row-major order と column-major order の違いの問題

(これはC言語に限った話ではないけれど…)

Fortran で

real a(3,2)

と宣言した場合、メモリーにおいては、 a(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2) の順に並ぶ (column-major order)。 MATLAB, Julia, Eigen (C++用のベクトル・行列用クラス・スライブラリィ) などは column-major order である。

ところが C で

double a[3][2];

と宣言した場合、メモリーにおいては、 a[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[][1], の順に並ぶ (row-major order)。 C++ や Python の NumPy などは row-major order である。

いわゆる線形演算をする場合、 LAPACK のような Fortran で書かれたサブルーチン・ライブライリィの利用を検討したくなる。その場合に、FORTRAN と C の2次元配列で、 order が異なっていることが問題となる。

同様に、2変数関数のグラフの鳥瞰図などを描く関数を呼ぶ場合も、それが問題となる。

個人的には、C や C++ でプログラムを書いていて、そこから Fortran で書かれたサブルーチンを呼び出す場合が多かったが、しばしば column-major で記憶するようにしたりした。

例えば column-major order の言語で書かれた関数 bird_view() を、 C で書いたプログラムから呼び出すために次のようなコード (2つ示す) を用いる。

方法1 を u[][] に記憶する

double u[ny+1][nx+1]; // ここの不自然な感じと ... for (i = 0; i <= nx; i++) { x = a + i * dx; for (j = 0; j <= ny; j++) { y = c + j * dy; u[j][i] = f(x, y); // ここの不自然な感じは目をつぶる } } bird_view(u, nx, ny);

方法2 マクロを使って column-major order でアクセス

// nx, ny がマクロ定義の場合でも正しく動くように (nx), (ny) と書く #define u(i,j) u[(i)+((nx)+1)*(j)] // column-major order int nx, ny; double *u; ... u = malloc(((nx)+1)*((ny)+1)); if (u == NULL) ... for (i = 0; i <= nx; i++) { x = a + i * dx; for (j = 0; j <= ny; j++) { y = c + j * dy; u(i,j) = f(x, y); } } bird_view(u, nx, ny);

それから、これは馬鹿馬鹿しいようだが、 bird_view() のような関数を呼ぶときに、転置した内容の変数を用意する、ということも考えられる。

小さなスタックをどうするかの問題

ところで件の講習会では、スタックに大きなデータを記憶するのが難しい、という問題を取り上げていた。例えば

#define Nx (1000) #define Ny (2000) int main(void) { double u[Nx][Ny]; // 1000*2000*8=16MB ... }

のようなプログラムがうまく動かない、ということなのだけど、これはクラシックな問題だよね。

今の macOS で普通にやるとスタックサイズがわずか 8MB ということだけれど、私がプログラミングを始めた頃は、 Cの某処理系のスタック・サイズのデフォールト・サイズは 2KB だったりしたのだ (完全に昔話モード)。

数値計算でその手の大きな変数は、 1つのプログラムでたくさん出てくるものではないので、 auto はやめて、static にするという手もある。

#define Nx (1000) #define Ny (2000) double u[Nx][Ny]; // 関数の外に出してしまう int main(void) { ... }

とか

#define Nx (1000) #define Ny (2000) int main(void) { static double u[Nx][Ny]; // static と指定する ... }

とか。

もっとも、この方法は、Nx, Ny が変数の場合には使えないので (可変長配列は static にできない)、あまり良いやり方ではない。

私なりの結論

結局は、以前からの結論を再確認することになりそう。

次が考えられる選択肢となる。

(a): malloc() を用いて1次元的に動的に確保する。ポインターのポインターを使って2次元配列風にアクセスする。
短所: (いつでも)少し余分なメモリーが必要になる。
(b): malloc() を用いて1次元的に動的に確保する。マクロを使って2次元配列風にアクセスする。
短所: 数が増えるとマクロを用意するのが面倒。
長所: 添字を 0 からでなくてから始めるとか、応用が利く。
(c): C言語以外のものを使う。C++ ならば C に近くて (そういう理由で選ぶべきではないと思うけれど)、 Eigen や Boost のような便利なクラス・ライブラリィがある。

(a)と(b)を比較すると、メモリー効率は (b) が優れているが、計算速度はケース・バイ・ケースのようである。

最近の私は (c) に傾きつつある。

個人的には、C はやはり設計の古い言語で、昔から言われて来たように高級なアセンブラに近いと思う。特にマクロとか使い出すと… (a) も (b) も理解できるけれど、人に強く勧める気にはなれない。 (アセンブラは重要で、必要性はなくなるはずがないと思うけれど、第一の選択肢にするものではない。)

C++ ならば、色々なクラス・ライブラリィが使いやすい。場合によっては、それが一番大きな理由となって、C++ の利用を勧めることがありうる。

桂田祐史
2020-04-20