GNU/Linux ユーザのための数値計算のTips

数値計算全般

C言語, C++, Fortran共通/混合

C言語, C++関係

Fortran関係

各種ツール,ライブラリ関係

LAPACK関係



計算結果が数値ではなく,NaNという文字列となってしまう.

このNaNとはNot a Numberの略語であり,0での除算を行った,あるいは,負の数のsqrtの計算を行ったなどの場合に発生します. (発生場所の特定方法)

計算結果が数値ではなく,Infという文字列となってしまう.

数値計算で通常用いられる倍精度の浮動小数点数(俗に倍精度実数)は,IEEE 754で定められているように64ビットのメモリに以下のように実数を格納しており,

先頭1ビット s 中間11ビット e 終端52ビット f
符号 指数部 小数部

メモリに格納された情報から実際の10進数での数値へは以下のように変換されます.

(-1)s × 2 e-1023 × F

ここで,Fは小数部fの52ビットを1.fのように並べて作った2進数を10進数に変換した実数です. ただし,指数部eのビットが全て0の場合,および全て1の場合は,

e f
全て 0 全て 0 0.0
少なくとも1つはゼロ以外 非正規数 (0.0と下限値の間を表現する)
全て 1 全て 0 Inf (無限大)
少なくとも1つはゼロ以外 NaN (非数)

であり,通常の数の表現では用いられません.

上記の表現形式のため,倍精度の浮動小数点数で表現できる値には,上限と下限があり,それぞれ以下となります.

上限 : 1.7976931348623157E+308
下限 : 2.2250738585072014E-308

この上限の値を超過すると数値計算上は無限大と扱われてInfとなります.

(上記の記述はSun Forte Developerのマニュアルの「数値計算ガイド」を参考にしました.) (発生場所の特定方法)

InfやNaNが発生する場所を特定したい.

演算の結果が上記のInfやNaNとなる時の検出のためにSigFPE(Floating Point Exception, 日本語では演算例外もしくは浮動小数点例外)というシグナルがあるのですが,数値計算を行う上では不合理なことに,デフォルトでは演算例外が起こってもSigFPEを発生せずに計算が続行されます.

そのため,InfやNanの発生場所の特定のためには,演算例外でこのSigFPEが発生するようにすればよいことになります. (すると,Kernel側でシグナルを検出してプログラムが停止します.) しかし,デバッグのためにプログラムをわざわざ書き換えるのは面倒なことです.

ところが,GNU glibcバージョン2.2以上とgcc(Gnu Compiler Collection)の組合せを使っているシステムでは,以下の手順によって,プログラムを一切書き換えることなく,演算例外でSigFPEを発生させてInfやNanの発生場所を特定することができます. (Windows+MinGWの場合は,こちらの手順になります.)

  1. 以下の内容のtrapfpe.cを記述する.
                #define _GNU_SOURCE 1
                #include <fenv.h>
    
                static void __attribute__ ((constructor))
                trapfpe ()
                {
                    feenableexcept (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
                }
              
  2. trapfpe.cを以下のようにコンパイルして,libtrapfpe.aを作成する.(root権限を持っているならば,さらにlibtrapfpe.aを/usr/local/libにインストールしておくと良い.)
                $ gcc -c trapfpe.c -o libtrapfpe.a
              
  3. NaNやInfの発生場所を特定したいプログラム(例えばsource1.cとsource2.c)をコンパイルするとき,以下のようにコンパイル時に"-g"オプションを指定し,リンク時に"-ltrapfpe"を指定して,デバッグオプション付きでlibtrapfpe.aをリンクした実行ファイルを作成する. (libtrapfpe.aを/usr/local/libにインストールしていない場合は,-Lオプションでlibtrapfpe.aがあるディレクトリを指定する必要がある. また,以下はC言語の例だが,C++やFortranではコンパイラにg++, g77, gfortranを用いて同様の手順を行えば良い.)
                $ gcc -c -g source1.c
                $ gcc -c -g source2.c
                $ gcc source1.o source2.o -lm -ltrapfpe
              
  4. 作成したa.outを引数としてgdbを実行する.
                $ gdb a.out
              
  5. (gdb) というプロンプトが出たら run と入力して,gdb内でプログラムを実行する. (実行時に引数が必要であれば,run の後ろに指定する.)
  6. InfやNaNが発生したらSigFPEが発生してプログラムが停止し,発生した場所の関数名,ソースファイル名と行番号,および発生した行が表示される.
    もし,これらの情報が表示されない時は,エラーを起こした数学ライブラリ(libm)の関数の中なので,(gdb) というプロンプトの後に up と入力して,その数学関数を呼んだ関数に制御を移すと,上記の情報が表示される.

なお,gdbはデバッガなので,変数の値を表示する("print 変数名")などデバッグに便利な機能をたくさん備えています. gdbの機能についてさらに詳しく知りたい人は,マニュアルなどを調べて下さい.

gdbを終了させるコマンドは"quit"です. "quit"を入力すると,

      The program is running.  Exit anyway? (y or n)
    

というメッセージが出力されて応答待ちとなるので,"y"を入力してgdbを終了させてください.

(上記のtrapfpe.cのプログラムはg77のinfoを参考にしました.)

また,gfortranの場合は以下のように実行ファイルを作成することでlibtrapfpe.aがなくても同じ効果が得られます. ただし,-ffpe-trapはリンクオプションではなくコンパイルオプションであることに注意してください.

      $ gfortran -c -g -ffpe-trap=invalid,zero,overflow source1.f
      $ gfortran -c -g -ffpe-trap=invalid,zero,overflow source2.f
      $ gfortran source1.o source2.o
    

配列の添字が定義(確保)した範囲から外れていないかどうか知りたい.

配列の添字が定義(確保)した範囲から外れることは,数値計算以外のプログラムにおいてバッファオーバーランと呼ばれ,UNIXにおけるセキュリティホールの主な原因となっております.

Fortranならばコンパイラのオプションで,配列の添字をチェックするコードの生成の有無が指定でき,コンパイラがg77でもgfortranでも,コンパイルオプションに"-fbounds-check"を指定するだけです. (添字のチェックはオーバーヘッドとなるので,デフォルトでは添字のチェックはしません.)

C言語やC++では,もともと数値計算のための言語ではなく,2次元配列のように見える配列の配列があるという事情のもと,多くのツールでは,各添字についてのチェックはせず,ポインタが配列として確保した領域を越えた外部にアクセスしようとした時に警告する という実装になっております. つまり,C言語での配列の添字の順序と合わせると,静的あるいは,Numerical Recipes in C方式で動的に多次元配列を連続した領域として確保した場合には,事実上1番目の添字のチェックのみとなります. 例えば,a[10][10]と確保した場合には,

a[-1][0] : NG
a[0][-1] : NG
a[-1][10] : OK (Numerical Recipes in C方式で多次元配列を確保した場合には NG となる)
a[0][10] : OK
a[1][-1] : OK
a[9][10] : NG
a[10][0] : NG

となり,a[-1][10], a[0][10], a[1][-1]では2番目の添字が定義した範囲を逸脱していますが,アドレスは確保した領域内になってしまうため,チェックできません. 逆にFortranは数値計算のための言語なので,上記のオプションでは,2次元配列(行列)でどちらの添字がはみ出た場合にも,エラーメッセージとともにプログラムが停止します.

C言語やC++で対象の配列が動的に確保されたものならば,コンパイラではなくライブラリのみで対応できるので色々なツールが開発されておりますが,GNU/Linuxではディストリビューションで用意されていることが多いElectric Fenceがシンプルで使いやすいのではないかと思います. このElectric Fenceの実体はlibefence.so/libefence.aであり,このライブラリをリンクしたプログラムでは,動的に確保された領域を越えて外部にアクセスしようとするとSegmentation Faultにより停止するので,コンパイル時に"-g"を指定し,リンク時に"-lefence"を指定して,上記のライブラリをリンクした実行ファイルを作り,gdbの中で動かすだけです. あとは,libtrapfpeの場合と同様です. なお,Electric Fenceのデフォルトでは,バッファオーバーラン,つまり配列の添字が定義した範囲を上回った場合のみしかチェックしません. 配列の添字が定義した範囲を下回る場合(バッファアンダーラン)についてチェックする場合には,環境変数EF_PROTECT_BELOWに値1を設定して実行します.

一方,C言語やC++で対象の配列が静的に確保されたものの場合には,コンパイラが対応する必要があります. gcc3.*までではそのような機能はありませんが,gcc4.0〜4.8ではMudflap(説明日本語訳)により,また,gcc4.8以降ではAddressSanitizerによって,動的および静的配列の領域チェック機能が提供されます.

Mudflapの使い方は,まずコンパイル時に"-g -fmudflap"を指定し,リンク時に"-lmudflap"を指定して,実行ファイルを作ります. この実行ファイルの挙動は環境変数MUDFLAP_OPTIONSにより変化するのですが,Mudflap自体のデフォルトでは配列領域を越えたアクセスに対してワーニングを出すだけでプログラムは停止しないので,環境変数MUDFLAP_OPTIONSの値に"-viol-segv"を指定し,領域を越えたアクセスでSegmentation Faultを発生させて,プログラムが停止するようにします. Mudflapのメッセージが過多と感じるならば,環境変数MUDFLAP_OPTIONSの値に"-viol-segv -no-verbose-violations"と設定して下さい. (環境変数MUDFLAP_OPTIONSに"-help"を指定して,実行するとオプションの詳細が得られます.) ただし,Mudflapはバッファアンダーランをチェックしないようです.

AddressSanitizerの使い方は,コンパイル時に"-g -fsanitize=address"を指定して,実行ファイルを作ります. この実行ファイルの挙動は環境変数ASAN_OPTIONSにより変化するのですが,AddressSanitizer自体のデフォルトでは配列領域を越えたアクセスに対してメッセージを出力してプログラムは正常終了するので,環境変数ASAN_OPTIONSの値に"abort_on_error=true"を指定して,領域を越えたアクセスでSigABRTによりプログラムが異常終了するようにします. AddressSanitizerのメッセージのうち,addressの出力を消すことは出来ないようですが,addressの出力の凡例は,環境変数ASAN_OPTIONSの値に"abort_on_error=true:print_legend=false"と設定すると消すことができます. (環境変数ASAN_OPTIONSに"help=true"を指定して,実行するとオプションの詳細が得られます.) なお,このAddressSanitizerは,バッファアンダーラン,バッファオーバーランの両方とも検出します.

そして,上記のようにMudflapあるいはAddressSanitizerを有効にしてコンパイルした実行ファイルをlibtrapfpeの場合と同様にgdbの中で動かせば,確保した領域を越えて外部にアクセスを試みた時点でプログラムがSigSEGVあるいはSigABRTにより停止します. ただし,停止した状態ではスタックの下部であり,カーネルやライブラリが提供する関数の内部にいるので,何回かupと入力して自分のソースファイルまで移動する必要があります. なお,Mudflapの場合では,自分のソースファイルまで戻った際には,エラーの発生したちょうどその行ではなく,少し場所がずれるようです. (スタックの性質上,関数がずれることは無いようですが.)

Fortranで記述されたライブラリをCのプログラムから利用したい.

Fortranの主な言語仕様は以下となります.

  1. Fortranでは返り値のない関数をサブルーチン(subroutine)と呼ぶ.
  2. Fortranプログラムには大文字と小文字は区別されない.
  3. Fortranプログラム内で HOGE() と定義されているサブルーチン/関数は,Cから見ると hoge_()のように,全て小文字となり最後にアンダースコアが付く.
  4. Fortranのサブルーチンと関数の引数は参照渡しである.

また,CとFortranの変数の型の対応の主なものは以下です.
(以下はlong intとintが同じ大きさとなる32ビットOSにおける代表的なものです. 64ビットOSに関しては以下と異なる場合がありますので,コンパイラのマニュアルで確かめて下さい. また,Fortranの3つある倍精度実数と倍精度複素数のうち,(KIND=8)がつくものはFortran90での宣言であり,g77では使うことができません.)

変数型 C Fortran
整数型 int INTEGER
倍精度実数型 double REAL*8
DOUBLE PRECISION
REAL (KIND=8)
文字型 unsigned char CHARACTER
論理型 int LOGICAL
倍精度複素数型 double _Complex
(C99 にて規定)
COMPLEX*16
DOUBLE COMPLEX
COMPLEX (KIND=8)

ここまでは,単純な置き換えだけで済む内容です. 実は,CとFortranの最大の差は,多次元配列,および文字列の取り扱いにあります. ただし,CとFortran間の文字列の受渡しは色々と面倒な内容を含むので,ここでは省略します. 必要な人は使用するFortranコンパイラのマニュアルで調べて下さい.

多次元配列についてのCとFortranの差は以下となります.

C Fortran
添字の下限値 0 1
データの並べ方 行優先(右側の添字が先に変化する) 列優先(左側の添字が先に変化する

例として2×2の行列 a を定義した場合を考えると,メモリ内には以下の順に配置されます.

C : a[0][0], a[0][1], a[1][0], a[1][1]
Fortran : a(1,1), a(2,1), a(1,2), a(2,2)

したがって,CからFortranのライブラリに行列を渡す場合は,データの転置が必要となります.

以上より,Cから整数,倍精度実数,整数,倍精度複素数,倍精度複素数行列の順で5つの引数をとるhoge5というFortranのサブルーチンを呼ぶ場合,C99 規格に準拠したコンパイラを用いるとすると,以下のようなプログラムとなります. (ただし,行列eの添字を自然な順序で用いるためには,転置操作が必要です.)

      void hoge5_(int *a, double *b, int *c, double _Complex *d, double _Complex *e);
      int a, c;
      double b;
      double _Complex d, e[10][10];

      hoge5_(&a, &b, &c, &d, &e[0][0]);
    

また,可読性は損なわれますが,このようなヘッダファイルを使い,10行10列の行列eの1行2列目の要素を参照するとき,e[idx2(1,2,10)]と書くことを許容するならば,以下のように転置操作がいらないプログラムを書くことができます. (なお,リンク先のヘッダファイルを読めば分かる通り,マクロ"FARRAY_DEBUG"を定義すると,上記のバッファオーバーランの検出法と組み合わせることを前提とした,バッファオーバーラン以外の添字チェックを行います.)

      #include "farray_index.h"

      void hoge5_(int *a, double *b, int *c, double _Complex *d, double _Complex *e);
      int a, c;
      double b;
      double _Complex d, e[10 * 10];

      hoge5_(&a, &b, &c, &d, e);
    

上記のプログラムをコンパイルするときは,Fortranライブラリを忘れずにリンクして下さい. (参照)

Fortranで記述されたライブラリをC++のプログラムから利用したい.

基本的には,上記のCからFortranのライブラリを呼ぶ場合と同じですが,C++にはCに加えて以下の特長があります.

  1. C++では仮引数にリファレンス型定数を指定できるため,呼ばれるライブラリ内で値を変更しないことが分かっている引数に対して変数を定義する必要はない.
  2. デフォルトで複素数型がある.
  3. クラスを用いることにより,Fortranライブラリに渡すときに転置を必要としない行列を定義できる.

上記の1. と2. の性質を使ったプログラムの具体例として,整数,倍精度実数,整数,倍精度複素数,倍精度複素数行列の順で5つの引数をとるhoge5というFortranのサブルーチンを呼ぶ場合で,最初の整数と倍精度実数の値がhoge5の中で変更されない場合は,以下のようなプログラムとなります. (ただし,行列eの添字を自然な順序で用いるためには,転置操作が必要です.)

      #include <complex>

      extern "C" void hoge5_(const int &a, const double &b, int *c,
      std::complex<double> *d, std::complex<double> *e);
      int c;
      std::complex<double> d, e[10][10];

      hoge5_(1, 2.0, &c, &d, &e[0][0]);
    

また,3. については,世の中にC++のために開発された行列クラスが沢山あり,中には行列の和や積が計算できるリッチなクラスもあります. ただ,一般的に言うと,リッチなクラスになればなるほど,ソースレベルで最適化をしなければクラスライブラリのオーバーヘッドが大きく,かといってソースレベルで最適化をすると最適化レベルの維持には環境(CPUやOS,コンパイラ,ライブラリなど)の更新に合わせたメンテナンスが必要となるという性質があります. そのため,リッチなクラスを使いたい場合には,メンテナンスがきちんと行われているものを選びましょう. (あるいは,自分できちんとメンテナンスをしましょう.)
それとは反対に,必要最低限のメソッドをinline関数でシンプルに実装するという アプローチのクラスもあります. こちらも色々な人が作っておりますが,例としてここに試作品を置いておきます. このクラスが何をしているか解読できるならば,使ってみて下さい. なお,このクラスでは行列eの1行2列目の要素にアクセスする際には,e(1,2)となります. このfmatrixクラスを用いると上記のhoge5を呼ぶプログラムは以下となります.

      #include <complex>
      #include "farray.hh"

      extern "C" void hoge5_(const int &a, const double &b, int *c,
      std::complex<double> *d, std::complex<double> *e);
      int c;
      std::complex<double> d;

      fmatrix<std::complex<double> > e(10, 10);
      hoge5_(1, 2.0, &c, &d, e());
    

なお,C++の規格ではtemplateのexportとしてinline関数ではない形で実装することも可能なはずなのですが,これができるC++コンパイラはあまりないので,templateを使うとinline関数で実装せざるを得ないという側面もあります.

インクルードするヘッダファイルの場所を-Iで指定するのが面倒.

自分で書いたヘッダファイルを各プログラムで共有するため"~/include"に置いた場合,gcc(Gnu Compiler Collection)は自動的に"~/include"の中を探すということはしないので,コンパイルオプションに"-I${HOME}/include"と指定する必要があります. しかし,以下の環境変数に対象のディレクトリを定義することにより,コンパイラの"-I"オプションで毎回指定するのと同じ効果が得られます.

環境変数 適用されるコンパイラ
CPATH gcc, g++
C_INCLUDE_PATH gccのみ
CPLUS_INCLUDE_PATH g++のみ

CやC++でπなどの各種定数を自分で定義するのが面倒.

数値計算プログラムをC言語とC++で書く場合,それぞれ<math.h>, <cmath>を必ずインクルードすると思いますが,glibcではこれらのヘッダファイル内に以下の各種定数が定義されています.

定数名
M_E e
M_LOG2E log2 e
M_LOG10E log10 e
M_LN2 loge 2
M_LN10 loge 10
M_PI π
M_PI_2 π/2
M_PI_4 π/4
M_1_PI 1/π
M_2_PI 2/π
M_2_SQRTPI 2/sqrt(π)
M_SQRT2 sqrt(2)
M_SQRT1_2 1/sqrt(2)

(ただし,gcc のオプションに "-std=c89" や "-std=c99" を指定すると,上記のマクロが使えないので,それぞれ "-std=gnu89" あるいは "-std=gnu99" とする必要があります.)

CやC++で大きな配列を定義すると,メモリは十分にあるはずなのに"Segmentation Fault"する.

CやC++では,普通に関数の内部(mainも含む)で変数を定義するとauto変数となり,stack領域から確保されます. しかし,通常ではこのstack領域の上限が設定されており,それ以上の記憶領域を確保しようとすると"Segmentation Fault"となります. (この上限はsh系のシェルならば"ulimit -a"で,csh系のシェルならば,"limit"で確認することができます.)

この"Segmentation Fault"を解決する方法には,以下の方法があります.

  1. 大きな配列を,"static"を指定した静的変数やグローバル変数として,stack領域以外の領域から確保させる.
    (この領域の上限は,sh系のシェルでは"ulimit -v"で表示されるvirtual memory, csh系のシェルでは"limit vmemoryuse"によりコントロールされる. 設定にもよるが,x86-linux の場合,通常はこの領域の上限はunlimitedになっている.)
  2. C言語ならばmallocやcalloc, C++ならばnewによって,配列を動的に確保する.
    (この場合にも,上限はulimit -v / limit vmemoryuseでコントロールされる.)
  3. stack領域の上限を変更する. sh系のシェルならば,
              $ ulimit -s unlimited
            
    csh系のシェルならば,
              # limit stacksize unlimited
            
    とすると,stack領域の上限を無効化できる. (ただし,あくまでも実装しているメモリの未使用部分しか実際に使うことはできませんが.)
    なお,上記のコマンドはシェルの内部コマンドであり,ログアウトするとデフォルトの制限に戻ります. 常にstack領域の上限を無効化したい場合には,以下の2つの方法があります.
    1. 個人レベルで解決する方法
      "~/.profile"や"~/.cshrc"に上記コマンドを記述する.
    2. ユーザ全体のレベルで解決する方法
      /etc/security/limits.confに以下の行を追加する.
                    *    soft   stack   unlimited
                  
      なお,古いPAM(Pluggable Authentication Modules)を使っているシステムではunlimitedが認識されず,新しくログインしても何も起こらないことがあります. 以下のコマンドで何も表示されないときは,そのシステムではunlimitedは認識されません.
                    # strings /lib/security/pam_limits.so | grep unlimited
                  
      そのような場合には,/etc/security/limits.confに上記の代わりに,以下のように指定するとカーネルの上限値で設定されます.
                    *    soft   stack   -1
                  
    このlimitの値にはsoft limitとhard limitの2種類があり,x86-linuxのstacksizeでは,通常hard limitがunlimitedでsoft limitのみが設定されているので,上記の記述をしております. (hard limitの値はそれぞれ,ulimit -Hs / limit -h stacksizeで確認できます.)
    使用しているシステムでhard limitも設定されている場合には,一般ユーザの権限ではhard limitの値を上げることができないので,/etc/security/limits.confに上記に加えてさらに以下も追加する必要があります.
              *    hard   stack   unlimited
            

なお,Windows+MinGWの場合は実行プログラム自体がstacksizeを設定しており,こちらの手順で変更します.

CやC++での真偽(0以外が真,0が偽という規則)で時々混乱する.

CやC++での真偽の規則は,シェルでのコマンドの戻り値と逆なので,無理に数字を覚えようとすると混乱します. コマンドの戻り値の方は変えられませんが,CやC++の方はbool型がありますので,真偽のみが必要な関数ならば戻り値にbool型を指定することで,数字を覚える必要が無くなります.

C++では"true"か"false"のどちらかの値を取る型として,bool型が標準規格に含まれています.

Cでのbool型はISO/IEC 9899:1999(通称 C99)にて規定された型であり,実体は_Bool型なのですが,ヘッダファイル<stdbool.h>にはbool型という別名と"true", "false"のマクロ定義が含まれており,C99に準拠したコンパイラならば,ヘッダファイル<stdbool.h>をインクルードするとC++と同様にbool型を使うことができます.

ただし,CやC++における暗黙の型変換により,bool型とint型の変数を比較する際にはbool型がint型に変換されるので,bool型のtrueはintの1に変換されて,下記のif文の中の...の部分は実行されません.

      int i = 2;
      if (i == true) {
          ...
      }
    

そのため,intを返す標準の組込み関数(例えばisalpha(), isdigit())により条件判断を行う場合には,下記のどれかとする必要があります.

  1. 下記のように,truefalseとの比較を行わず,if文の条件部に直接記述する.
              if (isalpha(c)) { ...
            
  2. 以下のように真偽のどちらの場合が必要なときも常にfalseとの比較を行い,trueとの比較を行わない.
              if (isalpha(c) != false) { ...
              if (isdigit(c) == false) { ...
            
  3. C言語ならば
              if ((bool)isalpha(c) == true) { ...
            
    C++ならば
              if (bool(isalpha(c)) == true) { ...
            
    のように明示的にキャスト変換を行う.

CやC++で真偽の2通りだけでなく,3通り以上の場合分けを作りたい.

C言語とC++には列挙という型がありますが,大抵の書籍では単に順番に定数を定義する方法として述べられていることが多く,利用価値を見出せない人が多いのではないかと思います. しかし,定数だけにとどまらず列挙型の変数や列挙型を返す関数が定義できるのが,一番大きなメリットではないかと筆者は考えます.

使いかたは下記の例のように,最初に列挙リストの名前と取り得る値のリストを定義して,変数や関数を定義するときにenumに加えて列挙リストの名前も指定します. 下記の例では使っていませんが,構造体のようにtypedefで新しい変数型のように宣言することもできます.

      enum trilogy {tubular_bells, hergest_ridge, ommadawn};
      enum trilogy favorite = ommadawn;
      enum trilogy todays_playlist(void);
    

列挙型の変数や関数を用いることによって,列挙のそれぞれがどのような値を取るかを気にする必要がなくなります. また,gcc/g++でコンパイルオプションに"-Wall"をつけた時,下記のようにswitchで列挙の全てを処理しておらず,さらに"default:"がないという条件が重なった場合に警告が出ます.

      switch (todays_playlist()) {
      case tubular_bells:
          make_free_minutes(50);
          break;
      case hergest_ridge:
          make_free_minutes(41);
          break;
      }
    

とはいえ,C言語では列挙型とint型の暗黙の型変換が相互に可能であり,列挙型変数に整数を代入してもエラーとならないので,気をつける必要があります. (C++では列挙型からint型への暗黙の型変換のみが可能なので,列挙型変数に整数を代入するとエラーとなります.)

C言語での複素数の取り扱いについて

現在,一般的に用いられているC言語は,最初にAmerican National Standard X3.159-1989として規定され,後にISO/IEC 9899:1990となった,通常"ANSI C"あるいは"C89"と呼ばれる規格のものです. この規格においては,数学関数は倍精度のみが定義され,複素数についてはまったく規定されていません. そのため,C言語による数値計算の入門書では,構造体で独自に複素数を実装していたりします. (独自にとは言っても,大体は似たような実装ですが.) しかし,C++のように演算子の多重定義ができないので,複素数を含んだ式がC++やFortranのような自然な式とならずに四則演算が関数呼び出しとなってしまい,数値計算におけるC言語の欠点と言われてきました.

その後ISO/IEC 9899:1999(通称 C99)によって,単精度とlong double(通常は4倍精度だが,倍精度と変わらないこともある)の数学関数と複素数型が規定され,最近のCコンパイラのほとんどはこのC99に準拠しています. (前述の_Bool型もC99で規定された型です.) このC99で,以下の複素数型とその四則演算,および整数型や実数型から複素数型への暗黙の型変換が規定されています.
("long int"があっても"int long"がないのと同じで,下記のように変数長の修飾子が先で,変数型の_Complexが後に来る方が正式なようなのですが,"_Complex double"などと逆にしても大丈夫なようです. また,純虚数を表す_Imaginary型も"optional"として規定されていますが,筆者にはこれをどのように使うのかがよく分かりません.)

この複素数型と共に,ヘッダファイル<complex.h>がC99で規定されており,この<complex.h>をインクルードすることで以下が使用可能になります.

_Complex の別名 complex
虚数単位 I, _Complex_I, _Imaginary_I(optional)
各種関数 csin, ccos, ctan, casin, cacos, catan, csinh, ccosh, ctanh, casinh, cacosh, catanh,
cexp, clog, cabs, cpow, csqrt, carg, cimag, creal, conj, cproj

上記の表の各種関数は,double _Complexに対するものだけを列挙しています. "conj", "cproj"以外は,複素数を示す頭の"c"を取ると,何をする関数かは自明と思います. "conj"は共役複素数を返す関数であり,"cproj"はリーマン球への射影を返す関数だということです.

<complex.h>内で"_Complex"の別名として"complex"が定義されていることより,以下の型が存在するのと同じことになります.
(先ほどのコメントと同様に,"complex double"などでも大丈夫なようです.)

また,実数から複素数への暗黙の型変換が規定されていることにより,以下のようなプログラムの記述が可能となります.

      #include <complex.h>

      double x = 1.0, y = 2.0;
      double complex z;

      z = x + I * y + clog(-10.0);
    

ただ,注意しなければならない点として,変数xが複素数の時に,例えばcexp(x)を書くつもりで間違えてexp(x)と書いた場合,xが暗黙の内に実数に変換されて,xの実数部のみの指数関数が計算されるという,非常に見付けにくいバグの原因となります. このため,保険として総称関数を有効にする<tgmath.h>もインクルードしておく方が良いでしょう. この<tgmath.h>をインクルードすると,上記の表の各種関数の欄のcarg, cimag, creal, conj, cproj以外の関数に対して,例えばexp(x)と書いたときに,変数xの種類によって実数版と複素数版を切替えるという,C++の関数のオーバーライドやFortranの総称関数と同じような役割を果たします. ただし,負の実数に対するsqrtlogで自動的に虚数を返すということはせず,<math.h>によるものと同じNaN(参照)を返すので,複素数が必要な時はcsqrtclogを使わなければなりません.
(当初は<tgmath.h>によるオーバーヘッドが大きいかと想像していましたが,色々と実験してみたところ,オーバーヘッドはかなり少ないようです. また,<tgmath.h>は自身で,<math.h>と<complex.h>をインクルードしているので,実はこれらの代わりに<tgmath.h>をインクルードするだけで十分です.)

一方,scanfprintfでの複素数の変換文字は用意されていないので,複素数を入力するときは実部と虚部を入れる実数型変数を用意して入力後に変換する必要があり,出力するときには,creal, cimagで実部,虚部を別々に実数型として出力しなければなりません.

C, C++, Fortranで複素数計算の簡単なベンチマークをしてみた所,C99はC++とほぼ同じ速さであり,libcの上にFortranライブラリのラッパがかかるFortranはやや遅いという印象を受けましたので,数値計算屋としてはC99がもっと世の中に普及して欲しいと思っております. しかし,JISやISOの規格書はありますが,一般向けにC99規格で書かれた書籍はまだまだ少ないのが現状なようです.
(追記:C++とC99の複素数演算について,以前はC++よりもC99の方が速いという認識でいましたが,良く調べてみるとC++はstreamによる出力が遅いのであり,複素数演算自体を比べるとC++とC99はほぼ互角でした.)

(C99の規格については,S. P. Harbison III and G. L. Steele Jr., C A Reference Manual, 5th edition, Prentice Hall, Upper Saddle River, New Jersey, 2002 を参考にしました.)

C言語での配列の動的確保について

書籍などで通常説明されているC言語であるISO/IEC 9899:1990(通称C89もしくはANSI C)では,配列を宣言する時その寸法は定数(#defineによるマクロ定数も含む)でなければなりませんでしたが,ISO/IEC 9899:1999(通称C99)において,可変長配列(variable-length arrays),つまり,const以外の変数でサイズを定義した配列の宣言が可能となりました. これは,今のところC99独自の規格であり,一般的なC++であるISO/IEC 14882:1998(E)(通称C++98)には無い機能です. (C++の新しい規格であるISO/IEC 14882:2003(E)(通称C++2003)ではどうかは調べていません.)

この可変長配列は,コンパイル時には寸法が不定であり,実行時に寸法が決まるという性質を持つため,malloc/callocと同様に実行時に動的確保されることとなります.

また,可変長配列と同時に,C++と同様の「任意の場所に変数宣言を記述可能」という性質もC99で追加されました.

これら2つのC99の新機能を組み合わせると,今までの以下のようなC89のプログラムの代わりに

      int N;
      double *array;

      scanf("%d", &N);
      array = (double *)calloc(N, sizeof(double));
    

C99では次のように書くことができます.

      int N;
      scanf("%d", &N);
      double array[N];
    

したがって,C言語初心者にとっての難関であるmalloc/callocを使わなくても良い場合が多いのではないかと思います. ただし,可変長配列はauto変数なので,デフォルトでは巨大な配列は宣言できないかも知れません. (参照)

また,可変長配列を使うとFortranでいう整合配列,つまり以下のような関数を記述することが可能になります.

      void flexible(int size, double array[][size]);
    

なお,C99の規格では引数をarray, sizeの順に並べることができないので,完全にFortran同等ではありません. しかし,gccでは以下の拡張構文により,flexible(array, size)と引数を並べる関数を記述することが可能です.

      void flexible(int size; double array[][size], int size);
    

上記のflexibleという関数の中で作業用配列が必要な場合には,やはり,可変長配列を用いて以下のように関数の中で普通に宣言できます.

      void flexible(int size, double array[][size])
      {
          double work[size][size];
          ....
      }
    

また,任意の場所に変数宣言を記述可能という性質について,今のところgccでは"-std=c99"あるいは"-std=gnu99"というオプションを必要としますが,C++と同様に

      for (int i = 0; i < 10; i++) {
          ....
      }
    

というfor文の記述が可能です.

(この項の記述は,以下を参考にしました.)

その他,C99の新機能について

bool変数, 複素数 可変長配列と任意の場所での変数宣言以外にC99で追加された新機能のうち,数値計算を行う上で役立ちそうなものとして以下があります.

  1. C++スタイルの"//"のコメント
  2. inline 関数
  3. double以外の浮動小数点数に対する数学関数
  4. <math.h>に追加された関数
  5. <tgmath.h>による総称関数

2.のinline関数はほぼC++と同じですが,C99では翻訳単位(普通は同一ファイル)で同じであれば良く,翻訳単位毎に違っていることを許容するので,結局 static 相当となるようです.

4.の<math.h>に追加された関数で,倍精度実数に関するものでは以下のようです.

丸め関連 nearbyint, rint, lrint, llrint, round ,lround, llround, trunc
剰余関連 remainder, remquo, ldexp, modf, scalbn, scalbln
指数・対数関数関連 exp2, expm1, log1p, log2, logb, ilogb
逆双曲線関数 asinh, acosh, atanh
特殊関数 erf, erfc, lgamma, tgamma
その他 cbrt, hypot, fma, fdim, fmax, fmin, signbit, copysign
IEEE 754 関連 isfinite, isinf, isnan, isnormal, nan, nextafter, nexttoward
条件判断 isgreater, isgreaterequal, isless, islessequal, islessgreater, isunordered

(この項の記述は,以下を参考にしました. また,http://gcc.gnu.org/c99status.htmlによると,gcc4ではinline関数,複素数,可変長配列のstatusが"Broken"だそうです.)

indentを用いた構文解析によるデバッグ

Fortranとは異なり,C言語やC++では言語仕様で関数や演算子の返す値が使われなくてもよいことになっています. (例えば,printf関数の返す値が使われることは,ほとんどありません.) この仕様のため,数式を入力する際に以下のように誤って"b"の後ろにセミコロンを入れても通常の構文として扱われ,"+ sin(c)"の部分は単に無視されてしまいます. このプログラムはgccのオプションに"-Wall"を指定してもワーニングとなりません. (sinのような関数呼び出しではなく単なる"+ c"ならば,"-Wall"オプションで検出されます.) また,GNU/Linuxでの標準的なlintであるsplintでは,下記のような簡単な例ならば検出できますが,自分の研究のために作成した数値解析プログラムでは検出できませんでした.

      a = b; + sin(c);
    

そこで,上記のような誤りを発見するために頼りになるのは,まずエディタのインデント機能です. エディタのインデント機能は行単位で行われるので,長い命令を複数行に分けた場合には意図しないインデントとして表れます.

あるいは,GNUのツールであるindentを用いると,上記のプログラムのような余分のセミコロンがある時は,indentにより予想外の出力結果が得られます. ただし,何もオプションを付けないでindentを動かすと,引数に指定したプログラムを書き換えてしまうので,自分のプログラムを書き換えたくない人はオプションに"-st"を指定して結果を標準出力に出力します. また,好みもあるとは思いますが,"-kr"オプションでKernighan & Ritchieスタイルに指定すると読みやすいのではないかと思います.

C++で大量のテキスト出力を行う場合の注意点

C++でプログラムを作るときは,せっかくだからC++らしくということで入出力にstream(<iostream>, <fstream>)を使いたくなると思いますが,このstreamは">>"と"<<"演算子を多重定義することによってどんなデータ型でも扱える汎用品なので,<cstdio>によるprintf()ファミリと比べると遅いです.
(三角関数の数表を作成するという簡単なプログラムでベンチマークを取ったところ,std::coutによる出力をprintf()による出力に変えただけでCPU時間が半分程度になりました.)

そのため,C++で計算結果の大量のテキスト出力がある場合にはC++らしさが薄れますが<cstdio>をインクルードしてprintf()ファミリを用いた方が,実行が速くなる可能性があります.
なお,入力についてもstreamとscanf()ファミリで同様となることが予想できますが,筆者はこちらの方は確認していません.

CとFortran9xの制御構文の比較

C Fortran
for (i = start; i <= end; i += incr) {
  ....
}
DO i = start, end, incr
  ....
END DO
while (condition) {
  ....
}
DO WHILE (condition)
  ....
END DO
do {
  ....
} while (condition)
N.A.
continue CYCLE
break EXIT
if (condition1) {
  ....
} else if (condition2) {
  ....
} else {
  ....
}
IF (condition1) THEN
  ....
ELSE IF (condition2) THEN
  ....
ELSE
  ....
END IF
switch (variable) {
case val1:
  ....
 break;
case val2:
  ....
 break;
default:
  ....
 break;
}
SELECT CASE (variable)
CASE (val1)
  ....
CASE (val2)
  ....
CASE DEFAULT
  ....
END SELECT
N.A. WHERE (logical_array)
  ....
ELSEWHERE
  ....
END WHERE
N.A. FORALL (i = start:end:incr)
  ....
END FORALL

Fortranで未定義の変数を使ったときに,エラーが出るようにしたい.

g77gfortranともにコンパイルオプションとして以下を指定するだけです.

コンパイラ コンパイルオプション
g77 -Wimplicit
gfortran -fimplicit-none

(サブルーチンと関数の各々に"implicit none"と指定しても同じですが, コンパイルオプションで対応した方が簡単です.)

Fortranで書いたプログラムの計算精度があまり良くない気がする.

C言語やC++とは異なり,Fortranでは "1.0" と書くと "default kind" つまり単精度の 実数定数として解釈されます. よって,通常は倍精度の計算が行われる現在では,精度を指定しない実数定数は計算精度の低下につながります. そこで,実数定数に精度を指定しなければならないのですが,これには2つの方法があります.

  1. "1.0d0" のように指数表示に"d"あるいは"D"を使う.
  2. "1.0_8" のようにアンダースコアの後に,"kind" を指定する.

上記の 1 は,Fortran77でも使用可能なのに対し,2 はFortran90以降のみとなります. この性質のため,例えば変数にπの値を与える場合には"pi = 4.0 * ATAN(1.0)" としてはならず,"pi = 4.0d0 * ATAN(1.0d0)"あるいは"pi = 4.0_8 * ATAN(1.0_8)"とする必要があります.

また,プログラム中で複素数を使っているならば,もう一つ気を付けなければならない点があります. 実数部 x, 虚数部 y の複素数を作る際に"CMPLX(x, y)" とやってしまうと,"default kind" の単精度で変換されます. 倍精度を指定して変換する方法は以下の2つとなります.

  1. "DCMPLX(x, y)" のように倍精度を指定する関数を使う.
  2. "CMPLX(x, y, kind = 8)" のように "kind" を指定する(Fortran9x以降のみ).

上記 1 の "DCMPLX" は正式な規格ではありませんが,Fortran77の拡張仕様としてほとんどの Fortran コンパイラで使用できます.

関数"REAL"は引数に複素数を与えて,複素数の実数部を取り出す関数として用いる場合のみ,返す値は同じ"kind"の実数となりますが,引数が整数や倍精度実数の場合には,返す値が単精度実数となりますので,混乱を避けるため,常に"kind"を指定するか,倍精度を指定する関数"DBLE"を使った方が良いでしょう.

一方,複素数の虚数部を取り出す関数 "AIMAG" は総称関数,つまり,引数の型によって自動的に関数の返す値の型が決まる関数であり,"kind" の指定は不必要で,指定するとエラーとなります.

なお,Fortran90/95/2003において,"DCMPLX"と"DOUBLE COMPLEX"は規格に入っておらず,"DOUBLE PRECISION"は推奨規格でないため,emacsのf90モードではキーワードとして認識されません. これらをキーワードとして認識させるためのlispファイルをここに置いておきます. このファイルを"~/.emacs.d"に置いた場合,"~/.emacs"に以下の行を追加することによって使うことが出来ます.

      (add-hook 'f90-mode-hook
                '(lambda ()
                   (load-file "~/.emacs.d/f90-keywords.el")))
    

Fortranの出力で改行を入れたくない.

Fortranでは,普通に "WRITE (6,*)" などで出力すると必ず最後に改行されます. Fortranの規格ではこの改行の抑制はできませんが,VAX/VMS 拡張 Fortran77ではFORMAT指定子の最後に "$" を指定すると,改行が抑制されます. この VAX/VMS 拡張 Fortran77 はいまは亡き DEC が VAX/VMS用のFortran77コンパイラのために設定した規格ですが,"do i = 1, 10 〜 end do"や"do while 〜 end do" など,ほとんどの Fortran77 コンパイラでサポートされており,事実上 Fortran77 の標準規格です. よって"$"Format指定子が,ほとんどのFortran77 コンパイラでサポートされているのはもちろんのこと,Fortran9x コンパイラはFortran77のプログラムをコンパイルできなければならないという事情のため,多くのFortran9xコンパイラでもサポートされています. (確か,最初に出た gfortran 4.0.0ではこのFormat指定子がサポートされていなかったような覚えがありますが,これを記述している時の最新の gfortran 4.0.2 ではサポートされています.)

これを使うと,

        WRITE (6, '(a,$)') 'Max number ? '
        READ (5, '(i5)') i
      

のように,入力を促すプロンプトを表示することができます.

Fortranプログラムでコマンド引数や環境変数を使用したい.

Fortranプログラムにおけるコマンド引数と環境変数の取り込みは,Fortran2003(ISO/IEC 1539-1:2004(E))でやっと規格に入りました.

コマンド引数の数は組込み関数"command_argument_count()"の返り値で知ることができ,引数そのものを知るには,組込みのサブルーチンを"call get_command_argument(i, arg)"のように使えば,文字変数"arg"がi番目の引数となります. また,"GALADRIEL"という環境変数の値を文字変数varに入れる場合には,"call get_environment_variable('GALADRIEL', var)"となります.

上記はあくまでもFortran2003の規格なので,gfortranでは引数に"-std=f2003"を付けないと"-Wall"を指定した場合に警告が出ます. また,Fortran2003に準拠していないコンパイラでは使うことができません. しかし,特にUNIXのコンパイラでは以前からこのような需要が多かったせいか,g77を含めた多くのFortranコンパイラでは拡張規格として,以下の関数/サブルーチンが用意されています.

Fortran 2003 拡張規格
command_argument_count() iargc()
get_commmand_argument(i, arg) getarg(i, arg)
get_environment_variable(name, value) getenv(name, value)

(Fortran2003の規格については,M. Metcalf, J. Reid and M. Cohen, Fortran 95/2003 explained, Oxford University Press, Oxford, 2004 を参考にしました.)

Fortranプログラムでπなどの各種定数を定義したい.

Fortran95以前ではparameterによる定数の初期設定において,sin()などの数学関数を使用することができませんでしたが,この機能もFortran2003においてようやく規格に採り入れられました.

そのため,Fortran2003では,以下の構文により各種定数を定義することができます.

      DOUBLE PRECISION, PARAMETER :: pi = 4.0d0 * ATAN(1.0d0)
      DOUBLE PRECISION, PARAMETER :: sqrt2 = SQRT(2.0d0)
    

Fortranプログラムで標準エラー出力(STDERR)を使いたい.

UNIXにおける標準I/Oストリーム入出力 STDIN, STDOUT, STDERR は通常,Fortranコンパイラにより以下のユニット番号に割り当てられています.

標準I/Oストリーム Fortranにおけるユニット番号
STDIN 5
STDOUT 6
STDERR 0

プログラムの実行に掛かった時間を計測したい.

UNIXではマルチユーザ/マルチタスクで動いておりCPUを 100% 使えるとは限らないので,他のプロセスの実行の状況によりプログラムの実行に掛かる実時間は大きく変化します. そのため,端末の前でプロンプトが返って来るまでの時間をストップウォッチなどで計測した結果は,実際の実行時間ではありません. また,実行に何時間も掛かるような大きなプログラムでは終了を待つのも大変です.

プログラム内の関数や実行文単位で実行時間を計測するためにはtimes(2)などのシステムコールを使うことになりますが,簡易的にプログラム全体の実行時間を計測すればよいのであればtimeコマンドがあります. timeコマンドの使いかたは簡単で,実行ファイルを引数としてtimeを起動するだけであり,実行ファイルが./a.outならば下記の例のようになります. (実行ファイルに引数が必要ならば,実行ファイルの後にそのまま指定するだけです.)

      $ time ./a.out
    

ただし,timeコマンド自体は大抵の場合シェルの組込みコマンドであるため,どのシェルを使っているかによって表示が異なります. (紛らわしいですが,外部コマンドとしてのtimeもあり,timeコマンドが組み込まれていない本当のBourne shellでは,sh系シェルの組込みコマンドとほぼ同じ動作である外部コマンドのtimeが使われます.) また,シェル組み込みのコマンドであることにより,下記のようにファイルにリダイレクトしたのでは./a.outの出力のみがoutput.txtに書き込まれてしまいます.

      $ time ./a.out > output.txt
    

したがって,計算時間をファイルにリダイレクトしたいのであれば,下記のようにサブシェル内で実行ファイルを動かす必要があります. (csh系のシェルの場合は後述するように標準出力に出力されるので>&ではなく,>でもよい.)

      $ (time ./a.out > output.txt) >& computation-time.txt
    

sh系のシェルを使っている場合には,実行プログラムの出力の後に,

      real    0m0.105s
      user    0m0.076s
      sys     0m0.024s
    

のように標準エラー出力に出力されます. 上記の各数字の意味は以下であり,

ラベル 意味
real プログラム実行の実際の時間(プロンプトが返って来るまでの時間)
user プログラムの実行に掛かったCPU時間
sys ファイルの読み書きなどでカーネル(OS)が消費したCPU時間

user, sysのCPU時間とはCPUの使用率を 100% として換算した時間です. 通常の場合,プログラムの実行時間として一番重要な時間はuserで示される時間となります.

csh系のシェルを使っている場合には,実行プログラムの出力の後に,

      0.076u 0.028s 0:00.10 90.0%     0+0k 0+0io 0pf+0w
    

のように標準出力に出力されます. 最初の2つの数字が先ほどのsh系のシェル組込みのtimeにおけるuser, sysに相当する時間(単位は秒)であり,最初の数字が一番重要な時間です. 3番目はsh系シェル組込みのtimeのrealに相当する時間ですが,これだけは"分:秒"のフォーマットです. (実行時間が1時間を越えた場合には,"時間:分:秒"のフォーマットとなります.) 4番目はCPUの使用率であり,(user + sys)/real で求めたパーセンテージです. 後ろの3つはカーネルが消費した実メモリ,ファイルの読み書きや仮想記憶に関する量なのですが,筆者は最初のuserに相当する時間のみしか気にしておりませんので,詳しいことは記述できません. 興味がある人は自分で調べてみて下さい.

ソースファイルが多くなったのでmakeを使ってコンパイルしたいが,Makefileを書くのが面倒,あるいはMakefileの文法が分からない.

世の中にはMakefileの文法を書いた本が色々と出版されていますので, Makefileを記述するのは難しいのではと後込みする人も多いと思います.

しかし,GNU makeを使うならば,デフォルトで以下のサフィックスルールが定義されており,Makefileを全く記述しなくてもファイルの拡張子から適切に判断して下表の右側のコマンドが実行されます. (複数の条件に当てはまるときは,上にある条件が優先されます.)

コマンド入力 条件 実行されるコマンド
make hoge.o hoge.c がある時 $(CC) $(CFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -c -o hoge.o hoge.c
hoge.cc がある時 $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -c -o hoge.o hoge.cc
hoge.C がある時 $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -c -o hoge.o hoge.C
hoge.cpp がある時 $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -c -o hoge.o hoge.cpp
hoge.f がある時 $(FC) $(FFLAGS) $(TARGET_ARCH) -c -o hoge.o hoge.f
hoge.F がある時 $(FC) $(FFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -c -o hoge.o hoge.F
make hoge hoge.o がある時 $(CC) $(LDFLAGS) $(TARGET_ARCH) hoge.o $(LOADLIBES) $(LDLIBS) -o hoge
hoge.c がある時 $(CC) $(CFLAGS) $(CPPFLAGS) $(LDFLAGS) $(TARGET_ARCH) hoge.c $(LOADLIBES) $(LDLIBS) -o hoge
hoge.cc がある時 $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(LDFLAGS) $(TARGET_ARCH) hoge.cc $(LOADLIBES) $(LDLIBS) -o hoge
hoge.C がある時 $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(LDFLAGS) $(TARGET_ARCH) hoge.C $(LOADLIBES) $(LDLIBS) -o hoge
hoge.cpp がある時 $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(LDFLAGS) $(TARGET_ARCH) hoge.cpp $(LOADLIBES) $(LDLIBS) -o hoge
hoge.f がある時 $(FC) $(FFLAGS) $(LDFLAGS) $(TARGET_ARCH) hoge.f $(LOADLIBES) $(LDLIBS) -o hoge
hoge.F がある時 $(FC) $(FFLAGS) $(CPPFLAGS) $(LDFLAGS) $(TARGET_ARCH) hoge.F $(LOADLIBES) $(LDLIBS) -o hoge

そのため,Makefileの文法として最低限必要なのは以下となります.

さらに,gccでは-MMオプションで,ソースファイル内の#include "hamatai.h"のように""で指定した#include文のみに対する依存関係の抽出ができます.
(正確には,-MMオプションは#include文で指定したヘッダファイルの依存関係のうち,/usr/include などのOS標準のディレクトリ内のファイルを除外するというものですが,通常のC言語の書き方に従うならば上記の動作になるはずです.)

以上を用いると,以下の手順でGNU make用の必要最低限のMakefileを記述することができます. なお,以下の手順によるMakefileを用いるためには,ソースファイルのうちのどれか1つのファイルの拡張子を取り除いたものを実行ファイルの名前とする必要があります.

  1. 自分が使っているプログラミング言語に合わせて下表のマクロで,コンパイラ,コンパイラのオプション,プリプロセッサのオプションを定義する.
    言語 コンパイラ コンパイルオプション プリプロセッサオプション
    C CC CFLAGS CPPFLAGS
    C++ CXX CXXFLAGS
    Fortran FC FFLAGS
  2. 実行ファイル作成時に必要なリンクオプションとライブラリをそれぞれ,LDFLAGS, LDLIBSに定義する. (コンパイルオプションの中でリンク時にも必要なものはLDFLAGSにも指定しておく.)
  3. オブジェクトファイルから実行ファイルを作成する時には常にCCが用いられるので,C言語以外のプログラミング言語でプログラムを記述しているならば,CCに自分の使っているコンパイラを指定する. (あるいは,LDLIBSに使用している言語のライブラリを指定する.)
  4. "実行ファイル:必要なオブジェクトファイル"という行を記述する. 例えば,実行ファイルorabidooを作成するためには,orabidoo.oに加えてfive.o, miles.o, out.oが必要ならば,
              orabidoo: orabidoo.o five.o miles.o out.o
            
    とする. 実際には"orabidoo.o"は指定してもしなくても同じ命令が実行されるのだが,陽に記述しておく方が分かりやすい. (実行ファイルを複数作成する場合には,ダミーターゲットを最初に指定する.例えば,tagoとmagoが作成したい実行ファイルの時は,
              all: tago mago
              tago: tago.o source1.o source2.o source3.o
              mago: mago.o source1.o source2.o source4.o
            
    とする.)
  5. エディタを一旦閉じて,ここまで作成した内容をファイルを書き出す. (最後の部分に,空行やコメントなどの目印を付けておくと良い.)
  6. 以下のようにソースファイルから依存関係を抽出して,Makefileに追加する. (以下は,Cの場合です. C++やプリプロセッサが必要なFortranの場合は,適宜変更して下さい. また,マクロを定義しなければならない場合は適宜追加する必要があります.)
              $ gcc -MM *.c >> Makefile
            

なお,上記の5, 6の手順の代わりにMakefileの末尾に下記を記述しても同じ動作となります. (ここでは,依存関係を指定するファイルを"make.dep"としています.)

      make.dep: *.c
      <TAB>gcc -MM $(CPPFLAGS) *.c > make.dep

      -include make.dep
    

ここで,<TAB>は,キーボードのタブキーの意味です.

以上の手順によるMakefileの典型例としてC版C++版を示します.

BSD由来のmakeの場合,TARGET_ARCHマクロとLOADLIBESマクロ,およびオブジェクトファイルから実行ファイルを作る暗黙のルールがなく,また,includeの対象ファイルを自動的に作成/更新する機能がありません. そのため,BSD由来のmakeを使う場合には,ヘッダファイルの依存関係の更新を自分で行い,実行ファイル作成の手順をMakefileに記述する必要があります.

また,gcc4がリリースされてまだ間もないせいもあるのだと思いますが,現状ではGNU makeにFortran9xのためのサフィックスルールが定義されておりません. そのため,gfortranでフリーフォームのFortran9x(拡張子: .f90/.F90)を使いたい人は,自分でサフィックスルールを記述する必要があります. 近いうちにGNU makeに取り込まれることは確実ですが,自分でサフィックスルールを記述するのが面倒な人はこちらをどうぞ. Makefileにインクルードするなり,MAKEFILES環境変数にフルパスで指定するなりすれば使うことができます. ただし,このルールはデフォルトのものよりも優先されますので,注意して下さい.

ソースファイルを色々な実行ファイルで共有しているのでライブラリを作りたいが,作り方がよく分からない.

現在のUNIXにおけるライブラリは,静的ライブラリ(static library)と共有ライブラリ(shared library)の2種類があり,それぞれ以下の特徴を持っています.

種類 ライブラリの作成 実行ファイルとの関係
静的 (static) 簡単 実行ファイルに組み込まれてしまうため,ライブラリを更新すると実行ファイルを作り直す必要がある
共有 (shared) バージョン番号があり,少し面倒 実行時にリンクされるため,実行ファイルはそのままでも更新したライブラリを使用できる

基本的には,オブジェクトファイル群から静的ライブラリ(libhoge.a)を作るためには,

      $ ar cru libhoge.a *.o
      $ ranlib libhoge.a
    

であり,共有ライブラリ(libhoge.so.2.0.1, libhoge.so.2, libhoge.so)を作るためには, オブジェクトファイルを作るときに-fPICオプション付きでコンパイルし,

      $ gcc -shared -Wl,-soname,libhoge.so.2 -o libhoge.so.2.0.1 *.o
      $ ln -s libhoge.so.2.0.1 libhoge.so.2
      $ ln -s libhoge.so.2 libhoge.so
    

という手順が必要です. しかし,この手順で静的と共有の両方のライブラリを作るためには,-fPICオプション付きでコンパイルしたオブジェクトファイルと,-fPICオプション無しでコンパイルしたオブジェクトファイルの両方を用意しなければなりません.

そのため,オープンソースのライブラリではオブジェクトファイルを入れるディレクトリを分けるなどの工夫をしているのですが,GNU libtoolを使うとより簡単に両方のライブラリを作成することができます. (もともと,libtoolは様々なUNIXでも共通のコマンドでライブラリを作成できるように開発されたラッパーなので,今回のように2種類のライブラリを作る場合に適用するのはまさに役不足なのですが...)

libtoolは,実はバージョンによってかなり動作が異なるのですが,バージョン1.5以上を用いてhoge1.cとhoge2.cからライブラリ libhoge.a, libhoge.so.2.0.1, libhoge.so.2, libhoge.so を作る場合には,

      $ libtool --tag=CC --mode=compile gcc -c hoge1.c
      $ libtool --tag=CC --mode=compile gcc -c hoge2.c
      $ libtool --mode=link gcc -rpath ${HOME}/lib -version-number 2:0:1 -o libhoge.la hoge1.lo hoge2.lo
      $ libtool --mode=install install libhoge.la ${HOME}/lib
    

とすると,共有オブジェクトとライブラリを .libs ディレクトリに分けながら,コンパイルし,ライブラリを作成して,${HOME}/libにインストールまで行います. 拡張子loとlaがそれぞれ,ダミーのオブジェクトとライブラリとなります. --tagにはそれほどの意味はないのですが,指定しなければエラーになるのと,現状ではCC, CXX, F77, GCJしか設定できないので,Cに対してはCC,C++に対してはCXX,Fortranに対しては一番妥当なF77を指定します. そして,--mode=compile以降に本当のコンパイラを指定します. ライブラリの作成の際に,インストールするディレクトリを指定しなければならないのが少し不思議なのですが,呪文だと思って下さい.

また,先ほどのMakefileの作成法を用いる場合には,コンパイラを指定するマクロとして以下を設定してオブジェクトファイルを作り,

言語 マクロ名
C CC libtool --tag=CC --mode=compile gcc
C++ CXX libtool --tag=CXX --mode=compile g++
Fortran FC libtool --tag=F77 --mode=compile g77 もしくは libtool --tag=F77 --mode=compile gfortran

最後に,

      $ libtool --mode=link gcc -rpath ${HOME}/lib -version-number 2:0:1 -o libhoge.la *.lo
      $ libtool --mode=install install libhoge.la ${HOME}/lib
    

と実行すればよいことになります.

自分でソースからコンパイルしたライブラリを/usr/local/libに置いたが,gccがそのライブラリを見付けることが出来ない.

システム標準のライブラリが入るディレクトリ(/lib, /usr/lib)以外のディレクトリにあるライブラリを使用するためには,/etc/ld.so.confでそのディレクトリを指定し,ldconfigでキャッシュファイル/etc/ld.so.cacheを更新する必要があります.

具体的には,root権限で /etc/ld.so.conf

      /usr/local/lib
    

という1行を加え,ldconfigを実行すればOKです.

また,root権限が得られず${HOME}/libにライブラリ(libhoge.so)をインストールした場合には,キャッシュファイルは必要なく,コンパイル時と実行時に${HOME}/libを探しにいく指定をすれば良いので,

      $ gcc hogehoge.c -Wl,-rpath,${HOME}/lib -L${HOME}/lib -lhoge
    

とコンパイルすればOKです.

連立1次方程式を解くためや固有値・固有ベクトルの計算などでLAPACKを使いたいのだが,どのサブルーチン/関数を使えばよいのか分からない.

LAPACKが提供する大量のサブルーチン/関数は,

  1. Driver Routines
  2. Computational Routines
  3. Auxiliary Routines

の3種類に分類されており,このうちで普通のユーザが使うことが想定されているのはDriver Routinesです. このDriver Routinesに対して,Quick Reference(postscript)や,検索プログラムが提供されているので,これらを使えば必要なサブルーチン名を特定することができると思います. (上記を見れば分かる通り,Driver RoutinesがさらにSimple DriversとExpert Driversに分かれていますが,普通はSimple Driversを使えば大丈夫です. なので,Quick Referenceの方では最初の2ページが重要な情報となります.) Quick Referenceでは,問題の分類に対して'S'で始まるサブルーチンと'C'で始まるサブルーチンの2つが記述されていますが,これらはそれぞれ単精度実数と単精度複素数のサブルーチン名であり,通常の倍精度の計算で使う際には頭の一文字をそれぞれ'D'と'Z'に変更します.

以上の手順で特定されたサブルーチンを使うためには,次に引数の意味を調べなくてはなりませんが,LAPACKがインストールされている場合には,マニュアルも通常はインストールされています. (ディストリビュージョンが提供するバイナリの場合には,マニュアルが別パッケージになっているかも知れません.) そのため,例えば一般非対称行列の固有値・固有ベクトルを倍精度で求める場合,上記の手順によりサブルーチン名が,実数計算の場合DGGEV,複素数計算の場合ZGGEVと特定されるので,あとは,UNIXのコマンドプロンプトで

      $ man dggev
    

あるいは

      $ man zggev
    

と入力すれば,各引数の解説を読むことができます. あとは,Fortranから使う場合には,このマニュアルに書かれているようにプログラムを記述すれば良く,CやC++の場合には,関数のプロトタイプ宣言と関数の呼び出しをそれぞれCの場合, C++の場合をもとにして記述すればよいことになります.

例えば,10行10列の複素非対称行列 A の固有値のみを求め配列 W に格納する場合の,CとC++のプログラム例は,それぞれ以下のようになります.

上記のプログラムでは,

としており,ldvl, ldvr は固有ベクトルを求めない場合でも1以上の値でないとエラーとなるため,1を指定しています.

LAPACKを使って逆行列を求めたい.

何故かLAPACKでは普通に逆行列を求めるサブルーチンがDriver Routinesの中に入っていないため,逆行列が必要なときは,

  1. 線形代数方程式を解くサブルーチンの右辺行列に単位行列を指定する.
  2. まず,LU分解やコレスキー分解などのサブルーチンを呼び出し,その次に,その結果を使って逆行列を計算するサブルーチンを呼び出す.

のどちらかの手順となります. メモリ使用量の面から見ると,右辺行列の領域を必要としない2.の手順の方が有利です.

上記の2.の手順に必要なサブルーチン名を特定するならば,検索アプレットの入口で"Computational Routines"の"Linear Equations"をクリックして出てくるアプレットで"Factor"と"Invert Using Factorization"を探して下さい.

行列計算でLAPACKを使うプログラムをC言語で書いたが,コンパイルすると"undefined reference to `e_wsfe'"のようなエラーメッセージが出て,実行ファイルを作成できない.

LAPACKはFortran77で書かれたライブラリであり,実行ファイル作成にはLAPACKをコンパイルしたFortranコンパイラのライブラリをリンクする必要があります. そのため,gcc2.*や3.*を使っている場合には,Fortran77を用いる予定がなくても,g77をインストールしてFortran77のライブラリを入れておきましょう. (gcc4.*の場合には,gfortranとなります. 他のコンパイラの場合は,対応するFortranコンパイラをインストールして下さい.)

g77でのFortranライブラリはlibg2c.aあるいはlibg2c.soであり,コンパイル時に-lg2cオプションを付けて,

      $ gcc hoge.c -llapack -lblas -lg2c -lm
    

のようにコンパイルすればOKです.

gcc4.*のFortranコンパイラgfortranでのFortranライブラリはlibgfortran.aあるいはlibgfortran.soとなり,コンパイル時のオプションは-lgfortranとなるので,

      $ gcc hoge.c -llapack -lblas -lgfortran -lm
    

のようにコンパイルすることになります.

なお,C言語から使うときはCLAPACKの方が良いと思う人もいるでしょうが,このCLAPACKNetlibのページに,"f2c'ed version of LAPACK"と書かれているようにf2cでC言語に変換したものであり,f2c方式がCとFortranのリンケージのデファクト・スタンダードとなった現在では,Fortranコンパイラを持っていなくてもLAPACKライブラリをビルドできることがその存在価値であり,C言語から使うためとしてはその役割を終えたと言ってよいでしょう. 実際,CLAPACKライブラリ自体の使いかたは,元のFortran77版のLAPACKと変わりません. そのため,gcc(Gnu Compiler Collection)を使っているならば, Fortranコンパイラは容易に入手できるので,CLAPACKを使う必要は全くありません. (特にdebianやRedHat系のディストリビュージョンでは,バイナリが提供されているので,自分でコンパイルする必要すらありません.) CLAPACKを自分でコンパイルする場合には,f2cが何をやっているかの知識が必要となることがあるので,個人的には,Fortran77版のLAPACKをコンパイルするよりも難しいのではないかと思っています.

一方,LAPACKに対するCLAPACKの利点として,ヘッダファイル"clapack.h"が付いてくるということがありますが,

  1. "clapack.h"を使うためには,"f2c.h"が必要.
  2. ISO/IEC 9899:1990(通称 ANSI C あるいは C89)で記述されているため,構造体として複素数を定義しなければならない.
  3. 大抵の場合は実際に呼び出すLAPACKの関数は非常に少ない.
  4. LAPACKの関数を使うためには,引数の意味や順序を調べることが必要.

ということより,自分が使うLAPACKの関数のプロトタイプ宣言をC, C++をもとにして自分で記述した方が色々と便利なのではないかと思います.

CPU に対して最適化されたBLAS,LAPACKを既存のプログラムで使いたいのだが,ソースがないので再コンパイル出来ない.あるいは,時間がかかるので再コンパイルしたくない.

例えば,Intel Math Kernel Libraryの場合ならば,BLASが入っているのがlibmkl.soであり,LAPACKが入っているのがlibmkl_lapack32.solibmkl_lapack64.soです. また,libmkl.soを使うためには,libguide.soもリンクする必要があります. 一方,コンパイル時に -llapack -lblas と指定したプログラムでは,

      $ ldd (対象のプログラム) | egrep '(blas|lapack)'
    

で表示されるライブラリ(現在のRedHat/Fedora Coreでは,libblas.so.3liblapack.so.3)を実行時に探します.

そのため,libblas.so.3liblapack.so.3という名前で,さらに上記の実体のライブラリを探しに行くファイルを作ればいいことになります.GNU ldではスクリプトが記述できるという特徴があるので,GNU/Linuxでは以下の手順となります. (以下では,対象のプログラムを octave とします.)

  1. 下記のように内容が1行のlibmkl_blas.so, libmkl_lapack.soというテキストファイルを作成する.
              $ echo "INPUT ( libmkl.so libguide.so )" > libmkl_blas.so
              $ echo "INPUT ( libmkl_lapack32.so libmkl_lapack64.so )" > libmkl_lapack.so
            
  2. 以下のようにlibblas.so.3, liblapack.so.3のバイナリを作成する.
              $ ld -shared -soname libblas.so.3 -o libblas.so.3 libmkl_blas.so
              $ ld -shared -soname liblapack.so.3 -o liblapack.so.3 libmkl_lapack.so
            
  3. 作成したlibblas.so.3, liblapack.so.3を適当なディレクトリ(例えば,/usr/local/mkl/lib)に入れる.
  4. 以下の内容のmkl-execというスクリプトを作成して,適当なディレクトリ(例えば,/usr/local/mkl/bin)に入れる.
              #! /bin/csh -f
              set mkl_utils_path = /usr/local/mkl
              setenv LD_LIBRARY_PATH ${mkl_utils_path}/lib
              set program = `where $0:t | egrep -v "(${mkl_utils_path}|\./)"`
              exec $program[1] $argv
            
  5. 作成したmkl-execスクリプトに実行属性を付ける.
  6. mkl-execスクリプトを入れたディレクトリ内で,
              $ ln -s mkl-exec octave
            
    とシンボリックリンクを張る.

上記の手順で作成されたoctaveは,もともとのBLAS, LAPACKではなく,Intel Math Kernel Libraryの提供するBLAS, LAPACKをリンクして実行します.

なお,Gnu Rでは,起動ファイル(Fedora Core Extraのrpmなら/usr/bin/R)がシェルスクリプトであり,その中で環境変数LD_LIBRARY_PATHを書き換えてしまうため,同様の手法を使うためにはRの起動ファイル中で定義されている変数R_LD_LIBRARY_PATHで,上記により作成したlibblas.so.3, liblapack.so.3が入っているディレクトリが,もともとのBLAS, LAPACKが入っているディレクトリよりも先になるように設定する必要があります.

また,もともとのBLAS, LAPACKを置き換えてしまいたい場合には,libblas.so.3, liblapack.so.3を上記のものに置き換え,さらに,libblas.so, liblapack.soをそれぞれ,上記のlibmkl_blas.so, libmkl_lapack.soで置き換えればよいのだと思いますが,筆者は確認しておりませんので,試してみたい方はリカバリーの手段を確保してから行ってみて下さい.

Last update: 2016.7.27

[機械力学研究室ホームページに戻る]   [数値計算に関する情報のページに戻る]