C, C++ による数値計算プログラミング
1998.10 初版
2007.11.7 最終改訂
私が普段使う数値計算の環境は PC-Unix 上で動作する fortran です。
しかし、今では PC-Unix を使う人は少数派であり、Windows で研究する人が多いでしょう。
Windows では fortran は非常にマイナーであるように思われます（cygwin 上で動作する f77
以外は高価であり、cygwin を使うくらいなら PC-Unix を使った方が楽）。
Windows では C や C++ のコンパイラがフリーで入手可能なので、
C や C++ で数値計算プログラムを組む人が多いようです。
Matlab で組む人もいるでしょう。
今では C, C++ に関してはコンパイラがフリーで手に入ります。
Linux などの Unix では C, C++ コンパイラがデフォルトで入っていますが、
Windows においても、Visual C++, Borland C++, cygwin 上の gcc などが
フリーのコンパイラとして入手可能です。
ここでは C, C++ で数値計算プログラムを組むことを考えます。
私の考えでは C は数値計算に向きませんが、C++ は非常に便利であり、
慣れれば fortran より開発効率が高いように思えます。
ただし、C++ はとても複雑な言語です。
まず、C が数値計算に向かない理由は以下の通りです。
- 複素数が扱えない
- 整合配列を扱うのが難しい。
上記の点は C++ を使うことにより解決することができます。それ
ぞれについて説明しましょう。
1. 複素数
私は電磁界理論の数値解析を行っていますが、電磁界では複素数を
使うことは日常茶飯事です。ところが、C では複素数型が使えませ
ん。struct を使って複素数を表す型を定義することは出来ますが、
加減乗除、sin, cos, exp などの演算が出来ません。加減乗除のた
めの関数を定義すると、いちいち足し算や掛算のたびに関数呼び出
しが発生し、速度が大幅に落ちると思われます。かといって
c = a * b という掛け算を
<![CDATA[
c.real = a.real * b.real - a.imag * b.imag
c.imag = a.real * b.imag + a.imag * b.real
]]>
といちいち書き下すとプログラムが非常に見にくくなります。
C++ を使うと上記の問題は解決できます。g++ 3.4.2 を使う場合
の例を示します。
<![CDATA[
#include <stdio.h>
#include <iostream>
#include <complex>
using namespace std;
int main()
{
complex<double> c1,c2,c3,c4;
double d1,d2,d3,d4;
c1 = complex<double>(1,2);
c2 = complex<double>(2,3);
c3 = c1 + c2;
cout << c1 << " " << c2 << " " << c3 << endl;
d1 = c1.real();
d2 = c1.imag();
d3 = abs(c1);
c4 = exp(c1);
cout << d1 << " " << d2 << " " << d3 << " " << c4 << endl;
}
]]>
2. 整合配列 ( 2 次元配列のサイズ指定 )
Fortran では
<![CDATA[
subroutine matrix(a,n)
dimension a(n,*)
]]>
のように配列のサイズをサブルーチンの引数として定義することが
出来ます。これを整合配列と言います。しかし、C では
<![CDATA[
sub(double a[][n])
]]>
のような宣言は出来ません ( n は引数 )。n を定数にすると( ex.
#define n 4 ) 上のような宣言も可能ですが、その場合 n が固定
値になってしまいます。例えば、連立一次方程式を解くサブルーチ
ンの場合、n はサブルーチンの呼び元で指定する必要があるので、
固定値では困ってしまいます。ゆえに、C では
<![CDATA[
sub(double a[],int n)
]]>
として、a_{ij} に対するアクセスは
<![CDATA[
a[i*n+j]
]]>
と書かなくてはなりません。これはかなり屈辱的です。
C 言語では a[i] は *(a+i) と等価であることを利用して、次のよ
うに書けば、Fortran のように a[i][j] の形式でアクセスできる
ようになります。
<![CDATA[
#include <stdio.h>
#include <stdlib.h>
main()
{
int i,j,k;
double a[3][4];
k = 1;
for(i=0; i<3; i++){
for(j=0; j<4; j++){
a[i][j] = k;
k++;
}
}
sub (a,3,4);
}
sub(double b[], int dim_i, int dim_j)
{
int i,j;
/* 以下の 3 行が必要 */
double **a;
a = (double **) malloc( sizeof(double*) * dim_i );
for(i=0; i<dim_i; i++) a[i] = b + i * dim_j;
for(i=0; i<dim_i; i++){
for(j=0; j<dim_j; j++){
printf("%8.2f ",a[i][j]);
}
printf("\n");
}
free(a); /* これも必要 */
}
]]>
上のようにすると a[i][j] の形式で配列にアクセスできますが、
配列 1 列分の余分のメモリ領域が必要であり、関数の先頭でアド
レスの代入が必要です。また、関数の終わりで free をする必要が
あります。これを忘れると、この関数を呼ぶたびにどんどんメモリ
を消費してしまいます。
これはかなり苦しいプログラムのように思われます。
Numerical Recipes in C ( ニューメリカルレシピ・イン・シー ) では
上記の関数側で行っていた作業を、配列を宣言する側で実行する方法を
採用しています。
<![CDATA[
/*
a[i][j] ---> *(a+i)[j] ---> *( *(a+i) + j )
address contents
-----------------
a b1
a+1 b2
a+2 b3
-----------------
b1 1
b1+1 2
b1+2 3
b1+3 4
b2 5
b2+1 6
b2+2 7
b2+3 8
b3 9
b3+1 10
b3+2 11
b3+3 12
*/
#include <stdio.h>
#include <stdlib.h>
main()
{
int i,j,k,dim_i,dim_j;
/* ここから配列作成 */
double **a,*b;
dim_i = 3;
dim_j = 4;
a = (double **) malloc( sizeof(double *) * dim_i );
b = (double *) malloc( sizeof(double ) * dim_i * dim_j );
for(i=0; i<dim_i; i++){
a[i] = b + i * dim_j;
}
/* 配列作成ここまで */
k = 1;
for(i=0; i<dim_i; i++){
for(j=0; j<dim_j; j++){
a[i][j] = k;
k++;
}
}
sub(a,3,4);
/* 配列使用後の後処理 */
free(b);
free(a);
}
sub(double **a, int dim_i, int dim_j)
{
int i,j;
for(i=0; i<dim_i; i++){
for(j=0; j<dim_j; j++){
printf("%8.2f ",a[i][j]);
}
printf("\n");
}
}
]]>
いずれにせよ、配列を使用する前に malloc して、使用後に free する
必要があるので、どう頑張っても C で整合配列を扱うのは難しそうに見えます。
なお、配列を宣言する側で malloc する場合、
<![CDATA[
for(i=0; i<dim_i; i++){
b = (double *) malloc( sizeof(double ) * dim_j );
a[i] = b;
}
]]>
とすると、データを入れる dim_i * dim_j 個の領域が、連続したアドレスに
確保することが保証されないように見えます。そうなると、
他の人が作った行列処理用ライブラリや fortran で作成したライブラリとの
接続が困難になるので、一気に dim_i * dim_j 個の
領域を確保するようにして下さい。
C++ を使うと整合配列を以下のように実現する
ことが出来ます。なお、この方法は本学科の 1998 年度生の松本く
んに教えて頂きました。この場を借りて厚く御礼を申し上げます。
<![CDATA[
#include <stdio.h>
#include <stdlib.h>
/*
[] をオペレータ定義 : a[i] = p_dat + i * dim_j
C 言語では : a[i] = *( a + i )
ゆえに、
a[i][j] ---> ( p_dat + i * dim_j )[j] ---> *( p_dat + i * dim_j + j )
~~~~ ~~~~~~~~~~~~~~~~~~~~~
*/
//---------------- 行列クラスの宣言 ------------
class matrix {
public:
matrix(int i, int j);
~matrix();
double *operator[](int i);
double &operator()(int i, int j);
int dim_i, dim_j;
int alloc_execute;
private:
matrix(const matrix&); // コピーコンストラクタをプライベートにして
double *p_dat; // 値渡しを禁止
};
// コンストラクタとデストラクタ
matrix::matrix(int i, int j)
{
p_dat = (double *) malloc( sizeof(double) * i * j );
if ( p_dat == NULL ){
fprintf(stderr,"Cannot alloc memory in matrix::matrix(i,j)\n");
exit(1);
}
alloc_execute = 1;
dim_i = i;
dim_j = j;
}
matrix::~matrix()
{
if ( alloc_execute == 1 ){
free(p_dat);
p_dat = NULL;
dim_i = 0;
dim_j = 0;
}
}
// a[i][j] 形式でアクセスする
inline double *matrix::operator[](int i) // インライン関数
{
return( p_dat + i * dim_j );
}
// a(i,j) 形式でアクセスする
inline double &matrix::operator()(int i, int j) // インライン関数
{
return( *(p_dat + i * dim_j + j) );
}
//-------------- サンプルプログラム ---------------
void sub(matrix &a) // 参照で受ける
{ // 値渡しにすると、コンパイルエラー
int i,j;
for( i = 0; i < a.dim_i; i++ ){
for( j = 0; j < a.dim_j; j++ ){
printf("%8.2f ",a[i][j]);
a(i,j) = a(i,j) * 2;
}
printf("\n");
}
}
main()
{
int i,j,k;
matrix a(3,4); // 静的に宣言
k = 1;
for( i = 0; i < a.dim_i; i++ ){
for( j = 0; j < a.dim_j; j++ ){
a[i][j] = k; // a(i,j) でも OK
k++;
}
}
sub(a);
for( i = 0; i < a.dim_i; i++ ){
for( j = 0; j < a.dim_j; j++ ){
printf("%8.2f ",a(i,j));
}
printf("\n");
}
}
]]>
matrix クラスでは inline 関数を定義しています。インライン関数
はヘッダファイルの中で定義して下さい。
ここで示したクラスには 2 つ問題があります。
- 値渡しで関数を呼び出すと、コピーコンストラクタが呼ばれ、
関数からリターンする時にデストラクタが呼ばれます。
c++ ではコピーコンストラクタを明示的に記述しない場合、コンパイラが
デフォルトのコピーコンストラクタを生成します。
デフォルトのコピーコンストラクタは変数をコピーするだけなので、
関数からリターンする時に呼ばれるデストラクタによって、
配列の領域が free されてしまいます。
この例では、コピーコンストラクタを private にすることで、
関数への値渡しを禁止しています。コピーコンストラクタの実装例は
後で示します。
- a = b のように、代入演算子を使うと、デフォルトの演算子 = では
変数をコピーするので、a と b の p_dat（データ領域）が同じ場所に
なってしまいます。従って、代入の瞬間に a.p_dat のアドレスが
消えてしまうので、その領域を free することが出来ません。
また、b.p_dat ( = a.p_dat ) はデストラクタが
呼ばれるたびに free されるので、2 回 free されることになり、
エラーが発生します。
ここでは double 型の配列を扱う例を示しましたが、
電磁界理論では複素数を扱うことが多いので、実際はテンプレート
を使って様々な型を取り扱えるようにする必要があります。
私が自分の研究で使っている行列クラスは これ です。
このクラスにはコピーコンストラクタと代入演算子も定義しています。
＜ 結論 ＞
以上の考察より、私には C で数値計算のプログラムを組むメリッ
トはあまり感じられません。むしろ、デメリットの方が多いように
思えます。
C++ は非常に便利なのですが、プログラミングの難易度は fortran よりは
遥かに高くなります。整合配列を使うだけでも、
「ポインタ」「コンストラクタ」「デストラクタ」「演算子のオーバーロード」などを
理解する必要があり、C++ の便利な機能を使うには「継承」「仮想関数」などの
理解も必要となります。これは、プログラミング初心者にはなかなか手ごわいです。
また、Fortran には信頼性の高い
数値計算ライブラリがありますが、C++ 用の数値計算ライブラリに
関しては未調査です。今では Matlab を使う人も多いので、
数値計算ライブラリなどというものは不要なのかも知れません。
なお、この文章は間違いを書いている可能性があります。間違いな
どがありましたら、お手数ですが、メールで御指摘下さい。