ソースコードを簡潔に直したいのですが…

質問

有限積分法(参考　http://www.akita-nct.jp/yamamoto/study/thesis/2005/thesis_namekawa.pdf)
を用いて、格子点 ( i, j, k )の磁場成分 Bx, By, Bz をその周辺を囲む電場成分 Ex, Ey, Ez で逐次計算するコードは次のようになります。
for ( i = 0; i < i_max; i++ ) {
　　for ( j = 0; j < j_max; j++ ) {
　　　　for ( k = 0; k < k_max; k++ ) {

　　　　　　Bx [ i ][ j ][ k ] = Ey [ i + 1 ][ j ][ k ] - Ey [ i + 1 ][ j ][ k + 1 ] + Ez [ i + 1 ][ j + 1 ][ k ] - Ez [ i + 1 ][ j ][ k ];

　　　　　　By [ i ][ j ][ k ] = Ez [ i ][ j + 1 ][ k ] - Ez [ i + 1 ][ j + 1 ][ k ] + Ex [ i ][ j + 1 ][ k + 1 ] - Ex [ i ][ j + 1 ][ k ];

　　　　　　Bz [ i ][ j ][ k ] = Ex [ i ][ j ][ k + 1 ] - Ex [ i ][ j + 1 ][ k + 1 ] + Ey [ i + 1 ][ j ][ k + 1 ] - Ey [ i ][ j ][ k + 1 ];

　　　　}
　　}
}

電場成分の係数行列 C ( 0, 1, -1 のいずれかをもつ ) を求める必要が生じた、すなわち

CCC　 E　 B
CCC * E = B
CCC　 E　 B

のような行列計算に変更するため、次のように書き直しましたが係数行列 C を求める部分を簡潔にできませんでした。

for ( i = 0; i < i_max; i++ ) {
　　for ( j = 0; j < j_max; j++ ) {
　　　　for ( k = 0; k < k_max; k++ ) {
　　　　　　for ( x = 0; x < 3; x++ ) {
　　　　　　　　row = index ( x, i, j, k, i_max, j_max, k_max );
　　　　　　　　y = ( x + 1 ) % 3;
　　　　　　　　z = ( x + 2 ) % 3;
　　　　　　　　if ( x == 0 ) {
　　　　　　　　　　C[ row ][ index ( y, i + 1, j,　　 k,　　 i_max, j_max, k_max )] =　1;
　　　　　　　　　　C[ row ][ index ( y, i + 1, j,　　 k + 1, i_max, j_max, k_max )] = -1;
　　　　　　　　　　C[ row ][ index ( z, i + 1, j + 1, k,　　 i_max, j_max, k_max )] =　1;
　　　　　　　　　　C[ row ][ index ( z, i + 1, j,　　 k,　　 i_max, j_max, k_max )] = -1;
　　　　　　　　}
　　　　　　　　if ( x == 1 ) {
　　　　　　　　　　C[ row ][ index ( y, i,　　 j + 1, k,　　 i_max, j_max, k_max )] =　1;
　　　　　　　　　　C[ row ][ index ( y, i + 1, j + 1, k,　　 i_max, j_max, k_max )] = -1;
　　　　　　　　　　C[ row ][ index ( z, i,　　 j + 1, k + 1, i_max, j_max, k_max )] =　1;
　　　　　　　　　　C[ row ][ index ( z, i,　　 j + 1, k,　　 i_max, j_max, k_max )] = -1;
　　　　　　　　}
　　　　　　　　if ( x == 2 ) {
　　　　　　　　　　C[ row ][ index ( y, i,　　 j,　　 k + 1, i_max, j_max, k_max )] =　1;
　　　　　　　　　　C[ row ][ index ( y, i,　　 j + 1, k + 1, i_max, j_max, k_max )] = -1;
　　　　　　　　　　C[ row ][ index ( z, i + 1, j,　　 k + 1, i_max, j_max, k_max )] =　1;
　　　　　　　　　　C[ row ][ index ( z, i,　　 j,　　 k + 1, i_max, j_max, k_max )] = -1;
　　　　　　　　}
　　　　　　}
　　　　}
　　}
}

int index ( const int d, const int i, const int j, const int k, const int i_max, const int j_max, const int k_max )
{
　　return ( ( d * i_max + i ) * j_max + j ) * k_max + k;
}

void Solver ( const int i_max, const int j_max, const int k_max)
{
　　int row, col,
　　　　cell_num = i_max * j_max * k_max;
　　for ( row = 0; row < cell_num; row++ ) {
　　　　for ( col = 0; col < cell_num; col++ ) {
　　　　　　B [ row ] += C [ row ][ col ] * E [ col ];
　　　　}
　　}
}

添え字計算が多く、非常に複雑なソースコードですが、規則性があるので、
皆様の力をお借りして簡潔に表現したいです。

よろしくお願いします。

回答No.3

Oh-Orange · ベストアンサー

★回答者 No.1 です。
＞x = 0のとき index関数の引数に i + 1、
＞x = 1のとき index関数の引数に j + 1、
＞x = 2のとき index関数の引数に k + 1
＞があるので、この部分をどうにかできればいいのですが…
　↑
　規則性に見えて規則的に処理できない気がします。多分。自信ないが…。
　それから for で繰り返しするよりは順番に記述したほうが早い気もします。
　前回の回答と同様に書き直すと以下のようになります。

サンプル:
int i0, j0, k0;
int i1, j1, k1;
int *pC;

for ( i1 = 1, i0 = 0 ; i0 < i_max ; i0++, i1++ ){
　for ( j1 = 1, j0 = 0 ; j0 < j_max ; j0++, j1++ ){
　　for ( k1 = 1, k0 = 0 ; k0 < k_max ; k0++, k1++ ){
　　　*pC = &C[ index(0,i0,j0,k0,i_max,j_max,k_max) ][ 0 ];
　　　*pC[ index( 1, i1, j0, k0, i_max, j_max, k_max ) ] = +1;
　　　*pC[ index( 1, i1, j0, k1, i_max, j_max, k_max ) ] = -1;
　　　*pC[ index( 2, i1, j1, k0, i_max, j_max, k_max ) ] = +1;
　　　*pC[ index( 2, i1, j0, k0, i_max, j_max, k_max ) ] = -1;

　　　*pC = &C[ index(1,i0,j0,k0,i_max,j_max,k_max) ][ 0 ];
　　　*pC[ index( 2, i0, j1, k0, i_max, j_max, k_max ) ] = +1;
　　　*pC[ index( 2, i1, j1, k0, i_max, j_max, k_max ) ] = -1;
　　　*pC[ index( 0, i0, j1, k1, i_max, j_max, k_max ) ] = +1;
　　　*pC[ index( 0, i0, j1, k0, i_max, j_max, k_max ) ] = -1;

　　　*pC = &C[ index(2,i0,j0,k0,i_max,j_max,k_max) ][ 0 ];
　　　*pC[ index( 0, i0, j0, k0, i_max, j_max, k_max ) ] = +1;
　　　*pC[ index( 0, i0, j1, k0, i_max, j_max, k_max ) ] = -1;
　　　*pC[ index( 1, i1, j0, k0, i_max, j_max, k_max ) ] = +1;
　　　*pC[ index( 1, i0, j0, k0, i_max, j_max, k_max ) ] = -1;
　　}
　}
}

int index( int d, int i, int j, int k, int i_max, int j_max, int k_max )
{
　return ((d * i_max + i) * j_max + j) * k_max + k;
}

void Solver( int i_max, int j_max, int k_max )
{
　int row, col, cell_num = i_max * j_max * k_max;
　int *pC = &C[0][0];

　for ( row = 0 ; row < cell_num ; row++ ){
　　for ( col = 0 ; col < cell_num ; col++ ){
　　　B[ row ] += *pC++ * E[ col ];
　　}
　}
}

その他:
＞y = ( x + 1 ) % 3;
＞z = ( x + 2 ) % 3;
　↑
　これや if で x の値を分岐するよりも順番に記述するだけでも良いと思います。
　ぱっと見て規則的に見えますが、これ以上の最適化は逆に見づらくなりませんか？
　index() 関数をマクロにして i_max、j_max、k_max の記述をなくすことも可能です。
　#define INDEX(d,i,j,k)　(((d * i_max + i) * j_max + j) * k_max + k)
　これなら
　*pC = &C[ INDEX(0,i0,j0,k0) ][ 0 ];
　*pC[ INDEX( 1, i1, j0, k0 ) ] = +1;
　*pC[ INDEX( 1, i1, j0, k1 ) ] = -1;
　*pC[ INDEX( 2, i1, j1, k0 ) ] = +1;
　*pC[ INDEX( 2, i1, j0, k0 ) ] = -1;
　とできます。単純に見やすいだけでも問題ない気がしますけど。
・以上。参考に。

お礼

再度ご回答ありがとうございます。

Oh-Orange様がご指摘のように
if 文で x の値を分岐させたのは苦肉の策でございまして…
確かに順番に記述した方が後々でも読みやすいですね。

index() 関数をマクロで書くのも参考になりました。

補足

あと気になった点がございまして、
ソートなどではデータを入れ替える代わりに
ポインタを張り替えるほうがよい…
とどこかで聞いたような気がします。

Oh-Orange様はポインタ演算で配列を参照していらっしゃいますが、
その方が高速になるなどのメリットがあるのでしょうか。
よろしくお願いいたします。

回答No.4

Oh-Orange

★返信。
・一般にポインタを使った方が配列の添え字からアドレスを計算しない分だけ
　早くなったりするのです。最近はコンパイラの最適化が優れているため
　row = index( 0, i0, j0, k0, i_max, j_max, k_max );
　C[ row ][ index( 1, i1, j0, k0, i_max, j_max, k_max ) ] = +1;
　C[ row ][ index( 1, i1, j0, k1, i_max, j_max, k_max ) ] = -1;
　C[ row ][ index( 2, i1, j1, k0, i_max, j_max, k_max ) ] = +1;
　C[ row ][ index( 2, i1, j0, k0, i_max, j_max, k_max ) ] = -1;
　としても機械語に変換されたときにはCPUのレジスタを使って *pC を
　行っているかもしれません。私は90年代からC言語を使っているため
　昔のコンパイラはあまり最適化が上手ではなく人間がポインタを積極的に
　使って高速化にしました。だからポインタを使いたがります。
　(去年とある質問から判明しましたがポインタを使うよりも添え字を使った方が)
　(実行処理が早かったことがあります。コンパイラの最適化が人間よりも優れて)
　(いるのだと実感しました。驚きです。ポインタの方が遅くなった。)
・以上。

お礼

>ポインタを使った方が配列の添え字からアドレスを計算しない
なるほど勉強になりました。

何度もご丁寧に回答して下さり、ありがとうございました。

回答No.2

jacta

詳しく見ていないので外しているかもしれませんが...
行列演算にBoost C++ LibrariesのuBLASを使えば、定義どおりの記述で済むのでは？

参考URL：

http://www.page.sannet.ne.jp/d_takahashi/boost/ublas/index.html

お礼

ご回答ありがとうございます。

質問の内容がうまく伝えられていないので誠に恐縮ですが、
単なる記述の簡潔さではなく、係数行列 C を求める部分の
( 0 で初期化した配列 C に -1, 1 のいずれかを代入する)
条件分けを簡潔に記述できたらと考えています。

補足

x = 0のとき index関数の引数に i + 1、
x = 1のとき index関数の引数に j + 1、
x = 2のとき index関数の引数に k + 1
があるので、この部分をどうにかできればいいのですが…

回答No.1

Oh-Orange

★一案として下さい。
・index()関数を使わずに次のようにしてみてはどでしょうか？
　単純に添え字カウンタを 0、+1 の２つを用意しています。
　後は Bx、By、Bz をポインタで操作しています。

char *pBx = &Bx[0][0][0];　←char型だと仮定(型がわからないので仮)
char *pBy = &By[0][0][0];
char *pBz = &Bz[0][0][0];
int i0, j0, k0;
int i1, j1, k1;

for ( i1 = 1, i0 = 0 ; i0 < i_max ; i0++, i1++ ){
　for ( j1 = 1, j0 = 0 ; j0 < j_max ; j0++, j1++ ){
　　for ( k1 = 1, k0 = 0 ; k0 < k_max ; k0++, k1++ ){
　　　*pBx++ = Ey[i1][j0][k0] - Ey[i1][j0][k1] + Ez[i1][j1][k0] - Ez[i1][j0][k0];
　　　*pBy++ = Ez[i0][j1][k0] - Ez[i1][j1][k0] + Ex[i0][j1][k1] - Ex[i0][j1][k0];
　　　*pBz++ = Ex[i0][j0][k1] - Ex[i0][j1][k1] + Ey[i1][j0][k1] - Ey[i0][j0][k1];
　　}
　}
}
以上。

お礼

ご回答ありがとうございます。

添え字カウンタとは妙案ですね。

Oh-Orange様は冒頭のソースに適用されていますが、
係数行列 C を求める部分にも使えそうです。

補足

ちなみに磁場成分 Bx, By, Bz (変更後は B ) と
電場成分 Ex, Ey, Ez (変更後は E )の型は
double 型で、係数行列 C は int 型です。