ポインタのポインタのポインタと3次元配列の利用

[2007.10.28]以下ははC言語での話。C++については01tec.Com - Online Payday Loans Direct Lenders Onlyをヒントにしてね。
[2007.6.16にかなり手を加えた]
流体の計算では、場を表現するために3次元配列を使う。ところがC言語で次の二つの条件を同時に満たすように書くのは難しい：

ある関数内、例えば最初の変数に関する差分を行う関数内では、配列の(i,j,k)番目の要素の呼び出しを a[i][j][k] と書きたい；
その関数では、どのようなサイズの3次元配列でも使えるようにしたい。配列のサイズを仮引数として与える形に書いておいて、独立したソースファイルとして(include あるいは分割コンパイルして)利用したい。

というのもプログラミング言語C 第2版 ANSI規格準拠 §5.7, p.137 にも書いてあるのだが

一般には、配列の最初の次元（添字）のみが自由で、他はすべて指定しなければならない

とある。同書に上げられている例としては、宣言

int daytab[2][13]

で作成された2次元配列に対して、これを関数 f に渡すときには f の宣言を

f(int daytab[2][13]) { ... }

f(int daytab[][13]) { ... }

f(int (*daytab)[13]) { ... }

のいずれかで書くとされている。これは配列の宣言、例えば int a[2][3][4]; が 2×3×4 個の連続したint型の配列要素を確保し、a がその先頭要素へのポインタでしかないからである。これでは汎用の関数を作成するのは難しい。
これに対する一つの解決法は大括弧 [ ] を用いた表記 "(なんとか)[i]" がポインタ(なんとか)で示される変数から i 番目の要素の内容、つまり "*((なんとか)+i)" であることを利用して、ポインタの連鎖使用による指示の列

a → a[i] → a[i][j] → a[i][j][k]

を作成することで、「配列要素を指示する」ようにすればよい。これで関数の仮引数の定義におけるサイズの明示的な指定を避けられる。この解決法では

表記 a[i][j][k] は *(*(*(a+i)+j)+k) を意味している。

これは配列 double a[n0][n1][n2] の a[i][j][k] が *(a+n0*n1*i+n1*j+k) であったのとは異なる。仕様とはいえ、なかなかに紛らわしいコードが書けるものである。
手許にあったC言語による数値計算入門―解法・アルゴリズム・プログラム (UNIX & Information Science)を参考に次のコードを書いたら、期待したとおりに動いた。ただしこのコードはメモリ確保に失敗したときの処理や、確保した領域の解放については書いていない。



#include <stdio.h>

#include <stdlib.h>

double ***malloc3d(int ix,int iy,int iz){

	double ***f;

	int i,j;

	f= (double***)malloc(ix*sizeof(double**));

	for(i=0;i<ix;i++){

		f[i]= (double**)malloc(iy*sizeof(double*));

	}

	for(i=0;i<ix;i++){ for(j=0;j<iy;j++){

		f[i][j]= (double*)malloc(iz*sizeof(double));

	}}

	return(f);

}

void sub0(double*** b){

	int i,j,k;

	for(i=0;i<3;i++){ for(j=0;j<3;j++){ for(k=0;k<3;k++){

		b[i][j][k]=100*i+10*j+k;

	}}}

}

void sub1(double*** b){

	int i,j,k;

	for(i=0;i<3;i++){ for(j=0;j<3;j++){ for(k=0;k<3;k++){

		printf("a[%d][%d][%d]=%f\n",i,j,k,b[i][j][k]);

	}}}

}

int main(void){

	double*** a= malloc3d(3,3,3);

	sub0(a);

	sub1(a);

	return(0);

}

結局、宣言文の double a[3][3][3] によるアロケーションとは似て非なる「配列もどき」である。見た目が３次元配列と同じなだけ。でも数値計算コードを書く身としては fortran とのアナロジーが簡単で便利。

2012.6.24 新コード

3次元データの保存領域が連続するようにアロケートする

#include <stdlib.h>
double *** mallocDouble3D ( int n0, int n1, int n2 ) {
  double ***b;
  int i,j;
  b= (double***) malloc( sizeof(double**) * n0 );
  for ( i= 0; i < n0; i++ ) {
    b[i]= (double**) malloc( sizeof(double*) * n1 );
  }
  b[0][0]= (double*) malloc( sizeof(double) * n0 * n1 * n2 );
  for ( i= 0; i < n0; i++ ) {
    for(j=0;j<n1;j++){
      b[i][j]= (double*) ( b[0][0] + ( j + i*n1 )*n2 );
    }
  }
  return b;
}
void freeDouble3D ( double ***b, int n0, int n1, int n2 ) {
  int i;
  free(b[0][0]);
  for ( i=0; i < n0; i++ ) {
    free( b[i] );
  }
  free(b);
}
int main (void) {
  double ***a= mallocDouble3D(2,3,4);
  // TODO's
  freeDouble3D(a,2,3,4);
}