コンパイルするファイル間の依存性はできるだけ減らそう 〜liboctaveへん〜

liboctave 利用時のコンパイルを「遅い」と感じたことはないだろうか．ここでは liboctave におけるヘッダファイルの依存関係を把握することにより，コンパイルを高速化する方法を検討する．
Scott Meyers (スコット・メイヤーズ): Effective C++ (吉川訳, アスキー出版局, 1998) が指摘しているように，コンパイルするファイル間の依存性はできるだけ減らした方がいい（34項）．特に liboctave のように大きなライブラリになると，コンパイル速度が顕著に違って来る． liboctave 利用時のコンパイルを高速化するために，ヘッダファイルの依存関係を把握しよう．

liboctave の場合， octave/config.h と octave/Matrix.h をインクルードしておけば，大概の機能を利用できるようになる． SVD (特異値分解) とか CHOL (コレスキー分解) さえも，利用できるようになる．これは裏を返せば，内部でファイルを大量にインクルードしているという事実に注意しなければならない．

例を示そう．

// cv1.cpp
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
using namespace std;

int main(int argc, char**argv)
{
  ColumnVector x(3,0.0), y(3,2.0);
  cout<<"x+y=\n"<<(x+y)<<endl;
  return 0;
}

このプログラムが依存しているヘッダを調べて欲しい． g++ だと，次のように -MM オプションをつければ，インクルードされているファイルを一覧できる：

g++ -MM -I/usr/include/octave-`octave-config -v` -L/usr/lib/octave-`octave-config -v` -loctave -lcruft cv1.cpp

この結果，次のようなリストが表示されるはずだ：

cv1.o: cv1.cpp /usr/include/octave-3.0.1/octave/config.h \
  /usr/include/octave-3.0.1/octave/oct-dlldefs.h \
  /usr/include/octave-3.0.1/octave/oct-types.h \
  /usr/include/octave-3.0.1/octave/Matrix.h \
  /usr/include/octave-3.0.1/octave/mx-base.h \
  /usr/include/octave-3.0.1/octave/MatrixType.h \
  /usr/include/octave-3.0.1/octave/boolMatrix.h \
  /usr/include/octave-3.0.1/octave/Array2.h \
  /usr/include/octave-3.0.1/octave/Array.h \
（長すぎるので中略）
  /usr/include/octave-3.0.1/octave/EIG.h \
  /usr/include/octave-3.0.1/octave/dbleLU.h \
  /usr/include/octave-3.0.1/octave/base-lu.h \
  /usr/include/octave-3.0.1/octave/CmplxLU.h \
  /usr/include/octave-3.0.1/octave/dbleQR.h \
  /usr/include/octave-3.0.1/octave/CmplxQR.h \
  /usr/include/octave-3.0.1/octave/dbleQRP.h \
  /usr/include/octave-3.0.1/octave/CmplxQRP.h

なんと，たった数行のプログラムのために，300を超えるヘッダがインクルードされていた．これは，コンパイルも遅くなるはずだ．

さて， ColumnVector は octave/dColVector.h で宣言されている．だから， octave/config.h と octave/dColVector.h さえインクルードしておけば， ColumnVector の基本的な機能は利用できるのである：

// cv2.cpp
#include <iostream>
#include <octave/config.h>
#include <octave/dColVector.h>  // 変更
using namespace std;

int main(int argc, char**argv)
{
  ColumnVector x(3,0.0), y(3,2.0);
  cout<<"x+y=\n"<<(x+y)<<endl;
  return 0;
}

もう一度， -MM オプションをつけて g++ を実行してみよう．今度は，

cv.o: cv2.cpp /usr/include/octave-3.0.1/octave/config.h \
  /usr/include/octave-3.0.1/octave/oct-dlldefs.h \
  /usr/include/octave-3.0.1/octave/oct-types.h \
  /usr/include/octave-3.0.1/octave/dColVector.h \
  /usr/include/octave-3.0.1/octave/MArray.h \
  /usr/include/octave-3.0.1/octave/Array.h \
  /usr/include/octave-3.0.1/octave/dim-vector.h \
  /usr/include/octave-3.0.1/octave/lo-error.h \
  /usr/include/octave-3.0.1/octave/lo-utils.h \
  /usr/include/octave-3.0.1/octave/oct-cmplx.h \
  /usr/include/octave-3.0.1/octave/syswait.h \
  /usr/include/octave-3.0.1/octave/MArray-defs.h \
  /usr/include/octave-3.0.1/octave/mx-defs.h

わずか10数個のヘッダしかインクルードされていないことがわかる．コンパイル速度を計測してみると，コンパイル時間に大きな差があることがわかる：

***% time g++ -I/usr/include/octave-`octave-config -v` -L/usr/lib/octave-`octave-config -v` -loctave -lcruft cv1.cpp
1.028u 0.080s 0:01.12 98.2%     0+0k 0+136io 0pf+0w
***% time g++ -I/usr/include/octave-`octave-config -v` -L/usr/lib/octave-`octave-config -v` -loctave -lcruft cv2.cpp
0.288u 0.048s 0:00.35 91.4%     0+0k 0+128io 0pf+0w

このように，後者だとコンパイル時間が1/3程度に短縮されていることが見て取れる．

使いたいクラスや関数がどのヘッダで宣言されているか調べるには， liboctave のソースを読めばよい．もしくは，ソースに doxygen を掛けて生成したドキュメント(HTML)をチェックすればするわかるはずだ．

以下，いくつかサンプルを示す．

#include <iostream>
#include <octave/config.h>
#include <octave/dColVector.h>
#include <octave/dRowVector.h>
#include <octave/dMatrix.h>
  // もともと octave/Matrix.h が include されていた
using namespace std;

int main(int argc, char**argv)
{
  Matrix A(2,3), B(2,2,1.0);
  ColumnVector x(3), y(2), z(2);
  for(int r(0); r<2; ++r)
    for(int c(0); c<3; ++c)
      A(r,c) = r+c;
  for(int r(0); r<3; ++r) x(r)=(r+1.0)*2.0;
  for(int r(0); r<2; ++r) y(r)=(r+1.0)*r;
  z = B*(A*x+y);

  cout<<"A="<<endl;
  cout<<A;
  cout<<"B="<<endl;
  cout<<B;
  cout<<"x= "<<x.transpose()<<endl;
  cout<<"y="<<y.transpose()<<endl;
  cout<<"z="<<z.transpose()<<endl;
  return 0;
}

これはかんたん．注意する必要があるのは， ColumnVector に transpose を適用すると RowVector になるから， octave/dRowVector.h もインクルードしておく必要があることくらいだ．

次の SVD (特異値分解) の例は，難しい．

#include <iostream>
#include <sstream>  // for stringstream
#include <octave/config.h>
#include <octave/dbleSVD.h>
#include <octave/mx-m-dm.h>
// #include <octave/mx-ops.h>

using namespace std;

int main(int argc, char**argv)
{
  Matrix A(3,6);
  stringstream ss (
    "  4   0   4   6   5   4"
    "  9   6   2   4   1   0"
    "  1   6   9   0   7   5"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // SVD オブジェクトの生成 :
  SVD  svd (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"svd.singular_values()= "<<endl<< svd.singular_values();
  cout<<"svd.left_singular_matrix()= "<<endl<< svd.left_singular_matrix();
  cout<<"svd.right_singular_matrix()= "<<endl<< svd.right_singular_matrix();
  // USV* を計算（Aに一致するか?） :
  cout<<"svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()= "<<endl
    << svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose();
  return 0;
}

SVD のために octave/dbleSVD.h をインクルードする必要があるのはすぐにわかるだろう．が，これだけだと，

svd.cpp:34: error: no match for 'operator*' in 'SVD::left_singular_matrix() const() * SVD::singular_values() const()'

というエラーが出る．つまり，

  Matrix operator* (const Matrix&, const DiagMatrix&);

が宣言されていないというわけ．これが宣言されているヘッダが octave/mx-m-dm.h なのだ．このヘッダを除いてみると，

// DO NOT EDIT -- generated by mk-ops
#if !defined (octave_mx_m_dm_h)
#define octave_mx_m_dm_h 1
#include "dMatrix.h"
#include "dDiagMatrix.h"
#include "mx-op-defs.h"
MDM_BIN_OP_DECLS (Matrix, Matrix, DiagMatrix, OCTAVE_API)
#endif

となっていて，さらに octave/mx-op-defs.h を見てみると，

#define BIN_OP_DECL(R, OP, X, Y, API) \
  extern API R OP (const X&, const Y&)
...
#define MDM_BIN_OP_DECLS(R, M, DM, API) \
  BIN_OP_DECL (R, operator +, M, DM, API); \
  BIN_OP_DECL (R, operator -, M, DM, API); \
  BIN_OP_DECL (R, operator *, M, DM, API);

と定義されていることから， octave/mx-m-dm.h をインクルードすると，

  extern OCTAVE_API Matrix operator + (const Matrix&, const DiagMatrix&);
  extern OCTAVE_API Matrix operator - (const Matrix&, const DiagMatrix&);
  extern OCTAVE_API Matrix operator * (const Matrix&, const DiagMatrix&);

という宣言が生成されることが分かる．人によるだろうが，結構ややこしく感じるだろう．

個別に operator を宣言しているヘッダを探すのが嫌な場合は， octave/mx-ops.h をインクルードすればいい．ただし，これによってコンパイル時間が倍以上になることを覚悟しなければならない．インクルードされるヘッダが一気に300近くに増加するからだ．残念ながら， octave/Matrix.h をインクルードした場合とほとんど同じ結果になってしまう．

逆行列はこんな感じでかんたん．

#include <iostream>
#include <octave/config.h>
#include <octave/dMatrix.h>
using namespace std;

void test1(void)
{
  Matrix A(3,3);
  A(0,0)=1.0; A(0,1)=2.0;  A(0,2)=3.0;
  A(1,0)=2.0; A(1,1)=-1.0; A(1,2)=1.0;
  A(2,0)=4.0; A(2,1)=3.0;  A(2,2)=2.0;
  int info;
  cout<<"A="<<endl<<A;
  cout<<"A.inverse(info)="<<endl<< A.inverse(info);
  cout<<"info= "<<info<<endl;
  cout<<"A.fill(0.0, 0,1, 2,1)="<<endl<< A.fill(0.0, 0,1, 2,1);
  cout<<"A.inverse(info)="<<endl<< A.inverse(info);
  cout<<"info= "<<info<<endl;
}

void test2(void)
{
  Matrix A(2,3);
  A(0,0)=1.0; A(0,1)=2.0;  A(0,2)=3.0;
  A(1,0)=2.0; A(1,1)=-1.0; A(1,2)=1.0;
  cout<<"A="<<endl<<A;
  cout<<"A.pseudo_inverse()="<<endl<< A.pseudo_inverse();
  cout<<"A * A.pseudo_inverse()="<<endl<< A * A.pseudo_inverse();
}

int main(int argc, char**argv)
{
  test1();
  cout<<endl;
  test2();
  return 0;
}

EIG (固有値分解) の例では， ColumnVector real (const ComplexColumnVector&) がどこにあるか探すのに手間取った．普通に octave/dColVector.h にあった．

#include <iostream>
#include <sstream>  // for stringstream
#include <octave/config.h>
#include <octave/EIG.h>
#include <octave/dColVector.h>  // for ColumnVector real (const ComplexColumnVector&);

using namespace std;

int main(int argc, char**argv)
{
  Matrix A(4,4);
  stringstream ss (
    "  4.0  2.0  3.0  1.0"
    "  0.0 -1.0  1.0  2.0"
    "  1.0  3.0  2.0  4.0"
    "  0.0 -1.0 -5.0 -4.0"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // EIG オブジェクトの生成 :
  EIG  eig (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // A を対称にする :
  A= 0.5*(A+A.transpose());
  cout<<"A="<<endl<<A;
  // 再度固有値分解 :
  eig= EIG(A,info);
  cout<<"info= "<<info<<endl;
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // 固有値，固有ベクトルを実数値で抽出する :
  cout<<"real(eig.eigenvalues()(2))= "<< real(eig.eigenvalues()(2))<<endl;
  cout<<"real(eig.eigenvectors().column(2))= "<<endl<< real(eig.eigenvectors().column(2));
  return 0;
}

よくつかうもの

私が普段よく利用するのは以下のライブラリだ：

#include <octave/config.h>
#include <octave/dColVector.h>
#include <octave/dRowVector.h>
#include <octave/dMatrix.h>
#include <octave/EIG.h>
#include <octave/dbleCHOL.h>
#include <octave/dbleSVD.h>
#include <octave/mx-m-dm.h>
// #include <octave/dbleDET.h>

これを全部インクルードしても， octave/Matrix.h をインクルードする場合に比べるとコンパイル時間が半分程度で済む（もっとも liboctave のヘッダのソースコード量がドミナントな場合の話だが）．