サイトの全体構成

liboctave のベクトル，行列，統計計算に関する機能

liboctave とは，Octave の機能（全機能ではないが）が集まった C++ 言語のライブラリファイルとヘッダファイル．自分で，C++ のプログラムを書き，liboctave の機能を呼び出すことが簡単にできる．

このページでは， ベクトル，行列，統計計算に関する機能のうちごく一部を説明する．

行列（基本操作、逆行列、固有値と固有ベクトル、SVD）

【お断り】以下，liboctave の機能を呼び出す C++ プログラムの見本と，コマンドラインでのビルドと実行の手順を説明する． 文法エラーは無いことは確認済みです．正しく計算できているかのチェックは済んでいません

前準備

Octave のビルドとインストールに関する Web ページの記述に従って，Octave のインストール
※ Octave 3.6 系列とそれ以前で liboctave の使用法に変化があったようで、このWebページでは、現時点での最新版である Octave 3.6 で説明を続けることにします．

(オプション) ソースコードのコピー

場合によってはソースコードがあった方が便利という場合があります（必ずしも必要というわけではありません）

cd /tmp
if [ ! -f octave-3.6.2.tar.gz ]; then 
    wget ftp://ftp.gnu.org/gnu/octave/octave-3.6.2.tar.gz
fi

sudo rm -rf octave-3.6.2
tar -xvzof octave-3.6.2.tar.gz
cd octave-3.6.2

cd /tmp/octave-3.6.2/src/DLD-FUNCTIONS
sudo cp *.cc /usr/local/include/octave-3.6.2/octave

コンパイルオブション

最初の見本は、行列のコンストラクタです．

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 3, 1.0);
  std::cout << X;

  return 0;
}

Ubuntu でのコンパイルオプションの例

上記のようなソースコードをa.cc のようなファイル名で保存しておく

ATLAS も Intel MKL も使わない場合
- コンパイルオプション: -O3 -fPIC -I/usr/local/include/octave-3.6.2 -I/usr/local/include
- リンクオプション: -O3 -fPIC -L/usr/local/lib/octave/3.6.2 -loctave -lcruft -loctinterp -llapack -lblas -lcblas -L/usr/local/lib -lcholmod -lumfpack -lcxsparse -lqrupdate -lfftw3 -lfftw3f -lreadline -lpcre -lgsl -lpthread -ldl -lm -lgfortran -lmetis -lamd -lcamd -lcolamd -lccolamd -lquadmath -lgomp
```
g++ -O3 -fPIC -I/usr/local/include/octave-3.6.2 -I/usr/local/include -c -o a.o a.cc
g++ -o a.out a.o -O3 -fPIC -L/usr/local/lib/octave/3.6.2 -loctave -lcruft -loctinterp -llapack -lblas -lcblas -L/usr/local/lib -lcholmod -lumfpack -lcxsparse -lqrupdate -lfftw3 -lfftw3f -lreadline -lpcre -lgsl -lpthread -ldl -lm -lgfortran -lmetis -lamd -lcamd -lcolamd -lccolamd -lquadmath -lgomp
```
ソースコードからビルドされた ATLAS と LAPACK を使う場合
- コンパイルオプション: -O3 -fPIC -I/usr/include/atlas -I/usr/local/include/octave-3.6.2 -I/usr/local/include
- リンクオプション: -O3 -fPIC -L/usr/local/lib/octave/3.6.2 -loctave -lcruft -loctinterp -L/usr/lib/atlas-base -ltmglib -llapack_atlas -llapack -lptcblas -lptf77blas -latlas -lcblas -L/usr/local/lib -lcholmod -lumfpack -lcxsparse -lqrupdate -lfftw3 -lfftw3f -lreadline -lpcre -lgsl -lpthread -ldl -lm -lgfortran -lmetis -lamd -lcamd -lcolamd -lccolamd -lquadmath -lgomp
```
g++ -O3 -fPIC -I/usr/include/atlas -I/usr/local/include/octave-3.6.2 -I/usr/local/include  -c -o a.o a.cc
g++ -o a.out a.o -O3 -fPIC -L/usr/local/lib/octave/3.6.2 -loctave -lcruft -loctinterp -L/usr/lib/atlas-base -ltmglib -llapack_atlas -llapack -lptcblas -lptf77blas -latlas -lcblas -L/usr/local/lib -lcholmod -lumfpack -lcxsparse -lqrupdate -lfftw3 -lfftw3f -lreadline -lpcre -lgsl -lpthread -ldl -lm -lgfortran -lmetis -lamd -lcamd -lcolamd -lccolamd -lquadmath -lgomp
```
64 ビット版の Intel Math Kernel Library を使う場合
- 環境変数: export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/intel/mkl/lib/intel64:/opt/intel/lib/intel64
- コンパイルオプション: -O3 -fPIC -m64 -DMKL_LP64 -I/opt/intel/mkl/include/intel64 -I/opt/intel/include/intel64 -I/opt/intel/mkl/include/fftw -I/usr/local/include/octave-3.6.2 -I/usr/local/include
- リンクオプション: -O3 -fPIC -m64 -DMKL_LP64 -L/opt/intel/lib/intel64 -L/opt/intel/mkl/lib/intel64 -L/usr/local/lib/octave/3.6.2 -loctave -lcruft -loctinterp -lmkl_intel_thread -lmkl_core -limf -lintlc -liomp5 -L/usr/local/lib -lcholmod -lumfpack -lcxsparse -lqrupdate -lfftw3 -lfftw3f -lreadline -lpcre -lgsl -lpthread -ldl -lm -lgfortran -lmetis -lamd -lcamd -lcolamd -lccolamd -lquadmath
```
g++ -O3 -fPIC -m64 -DMKL_LP64 -I/opt/intel/mkl/include/intel64 -I/opt/intel/include/intel64 -I/opt/intel/mkl/include/fftw -I/usr/local/include/octave-3.6.2 -I/usr/local/include  -c -o a.o a.cc
g++ -o a.out a.o -O3 -fPIC -m64 -DMKL_LP64 -L/opt/intel/lib/intel64 -L/opt/intel/mkl/lib/intel64 -L/usr/local/lib/octave/3.6.2 -loctave -lcruft -loctinterp -lmkl_intel_thread  -lmkl_core -limf -lintlc -liomp5 -L/usr/local/lib -lcholmod -lumfpack -lcxsparse -lqrupdate -lfftw3 -lfftw3f -lreadline -lpcre -lgsl -lpthread -ldl -lm -lgfortran -lmetis -lamd -lcamd -lcolamd -lccolamd -lquadmath 
```

行列（基本操作、行列と行列の積、LU分解、逆行列、行列式、SVD、QR分解、固有値と固有ベクトル、など）

コンストラクタ

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 3, 1.0);
  std::cout << X;

  return 0;
}

要素の参照

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 3, 1.0);
  X(0,0) = 1; X(0,1) = 3; X(0,2) = 5; 
  X(1,0) = 2; X(1,1) = 4; X(1,2) = 6; 
  std::cout << X;

  return 0;
}

行ベクトル，列ベクトルの取得

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 3, 1.0);
  X(0,0) = 1; X(0,1) = 3; X(0,2) = 5;
  X(1,0) = 2; X(1,1) = 4; X(1,2) = 6;
  
  RowVector r0 = X.row(0);
  RowVector r1 = X.row(1);
  ColumnVector c0 = X.column(0);
  ColumnVector c1 = X.column(1);
  ColumnVector c2 = X.column(2);
  
  std::cout << "r0=" << r0 << "\n";
  std::cout << "r1=" << r1 << "\n";
  std::cout << "c0=" << c0 << "\n";
  std::cout << "c1=" << c1 << "\n";
  std::cout << "c2=" << c2 << "\n";
  return 0;
}

行数，列数

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 3, 1.0);
  
  int nrows = X.rows();
  int ncols = X.cols();
  
  std::cout << X;
  std::cout << "nrows = " << nrows << "\n";
  std::cout << "ncols = " << ncols << "\n";
  return 0;
}

行列と行列の積

※ dMatrix.h

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;

  Matrix Y(2, 2, 1.0);
  Y(0,0) = 1; Y(0,1) = 2;
  Y(1,0) = 3; Y(1,1) = 4;

  Matrix Z = X * Y;

  std::cout << X;
  std::cout << Y;
  std::cout << Z;

  return 0;
}

LU 分解 (LU decomposition)

※ dbleLU.h

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;

  LU lu(X);
  Matrix L = lu.L();
  Matrix U = lu.U();
  PermMatrix P = lu.P();

  std::cout << L;
  std::cout << U;
  std::cout << Matrix( P );

  return 0;
}

逆行列

※ dMatrix.h

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;

  
  Matrix Z = X.inverse();
  
  std::cout << Z;
  std::cout << Z.inverse();

  return 0;
}

行列式 (determinant)

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;
  
  DET det = X.determinant();
  
  std::cout << X;
  std::cout << det.value() << "\n";

  return 0;
}

正方行列の Singular Value Decomposition (SVD) (singular value decomposition of A, A = U * S * D)

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;

  
  SVD svd(X);
  Matrix     U = svd.right_singular_matrix();
  DiagMatrix S = svd.singular_values();
  Matrix     V = svd.left_singular_matrix();
  
  std::cout << U;
  std::cout << S;
  std::cout << V;
  std::cout <<  (U * S * V).transpose();

  return 0;
}

QR 分解 (QR factorization)

※ dbleQR.h

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;

  QR qr(X); 

  std::cout << X;
  std::cout << qr.Q() << "\n";
  std::cout << qr.R() << "\n";
  std::cout << qr.regular() << "\n";  

  return 0;
}

正方行列の固有値と固有ベクトル (eigenvalues and eigen vectors of a matris)

#include<iostream>
#include<stdio.h>
#include<time.h>
#include<octave/config.h>
#include<octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 2, 1.0);
  X(0,0) = 3; X(0,1) = 5;
  X(1,0) = 7; X(1,1) = 9;

  
  EIG eig(X, /* calc_eigenvectors */ true );
  ComplexColumnVector v = eig.eigenvalues();
  ComplexMatrix L = eig.eigenvectors();
  
  std::cout << v;
  std::cout << L;
  std::cout <<  L * ComplexMatrix(ComplexDiagMatrix(v)) * L.inverse();

  return 0;
}

Octave の DLD-FUNCTIONS を C++ から使いたい場合

使いたいライブラリファイルを探す
ここの説明では，「/usr/local/lib/octave/3.6.2/rand.a」を使いたいとする
使いたい関数を探す
◆ 操作手順の例
```
nm -C /usr/local/lib/octave/3.6.2/rand.a 
```
◆ 実行結果の例

ここの説明では，「Frand(octave_value_list const&, int)」を使いたいとする。

ソースコードを書く

要点は３つ

「extern octave_value_list Frand(octave_value_list const&, int)」を入れる
DLD-FUNCTIONS の関数の引数については、octave_value_list を使うのが便利．
関数の返り値は octave_value_list 型の変数（例えば result）に入れ、result(0) のようにして取り出す

◆ ソースコードの例

// liboctave を用いた 乱数
#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/defun-dld.h>

#define DIM 3

extern octave_value_list Frand(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(DIM, DIM);
  std::cout << X;

  octave_value_list L;
  L.append(DIM);
  L.append(DIM);
  octave_value_list result = Frand(L, L.length());
  X = result(0).float_array_value (); 
  std::cout << X; 

  return 0;
}

ビルド
＜要点＞
リンクオプションとして /usr/local/lib/octave/3.6.2/rand.a を付ける
※ コンパイルオプションについては、この Web ページの上のほうで説明している．
◆ 実行結果の例

動作確認

./a.out

◆ 実行結果の例

コンボリューション (DLD_FUNCTIONS の Fconv2)

conv2.a の中の Fconv2() を使う。

◆ ソースコードの例

// liboctave を用いた 2次元畳み込み（コンボリューション）
#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/defun-dld.h>

#define DIM 3
#define DIM2 2

extern octave_value_list Fconv2(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(DIM, DIM);
  Matrix B(DIM2, DIM2);
  Matrix Z(DIM, DIM);

  X(0,0) = 1; X(0,1) = 4; X(0,2) = 7; 
  X(1,0) = 2; X(1,1) = 5; X(1,2) = 8; 
  X(2,0) = 3; X(2,1) = 6; X(2,2) = 9; 
  B(0,0) = 1; B(0,1) = 1; 
  B(1,0) = 1; B(1,1) = 1; 
  std::cout << X;
  std::cout << B;
  fprintf( stderr, "Z = Fconv2(X, B, \"full\")\n" );


  octave_value_list L;
  L.append(X);
  L.append(B);
  L.append("full");
  octave_value_list result = Fconv2(L, L.length());
  Z = result(0).float_array_value (); 
  std::cout << Z; 

  return 0;
}

◆ 実行結果の例

一様乱数 (DLD_FUNCTIONS の Frand)

rand.a の中の Frand() を使う。

◆ ソースコードの例

// liboctave を用いた 乱数
#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/defun-dld.h>

#define H 8
#define W 3


extern octave_value_list Frand(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(H, W);
  std::cout << X;

  octave_value_list L;
  L.append(H);
  L.append(W);
  octave_value_list result = Frand(L, L.length());
  X = result(0).float_array_value (); 
  std::cout << X; 

  return 0;
}

◆ 実行結果の例

正規乱数 (DLD_FUNCTIONS の Frandn)

分布が正規分布になるような乱数. rand.a の中の Frandn() を使う。

◆ ソースコードの例

// liboctave を用いた 乱数
#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/defun-dld.h>

#define H 8
#define W 3


extern octave_value_list Frand(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(H, W);
  std::cout << X;

  octave_value_list L;
  L.append(H);
  L.append(W);
  octave_value_list result = Frandn(L, L.length());
  X = result(0).float_array_value (); 
  std::cout << X; 

  return 0;
}

◆ 実行結果の例

1次元のFFT (DLD_FUNCTIONS の fft)

fft.a の中の Ffft() を使う。

◆ ソースコードの例

#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/ov-flt-complex.h>
#include<octave/ov-complex.h>
#include<octave/fCMatrix.h>
#include<octave/CMatrix.h>
#include<octave/defun-dld.h>

#define LEN 256

extern octave_value_list Frandn(octave_value_list const&, int); 
extern octave_value_list Ffft(octave_value_list const&, int); 
extern octave_value_list Fifft(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  ComplexMatrix X(1, LEN);
  ComplexMatrix Y(1, LEN);

  octave_value_list L;
  octave_value_list result;

  L.append(1);
  L.append(LEN);
  result = Frandn(L, L.length());
  X = result(0).complex_array_value ();  
  std::cout << "X\n"; 
  std::cout << X; 

  L.clear();
  L.append(X);
  L.append(LEN);
  result = Ffft(L, L.length());
  Y = result(0).complex_array_value ();  
  std::cout << "Y\n"; 
  std::cout << Y; 

  return 0;
}

◆ 実行結果の例

2次元のFFT (DLD_FUNCTIONS の fft)

fft.a の中の Ffft2() を使う。

◆ ソースコードの例

#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/ov-flt-complex.h>
#include<octave/ov-complex.h>
#include<octave/fCMatrix.h>
#include<octave/CMatrix.h>
#include<octave/defun-dld.h>

#define SZ 16

extern octave_value_list Frandn(octave_value_list const&, int); 
extern octave_value_list Ffft2(octave_value_list const&, int); 
extern octave_value_list Fifft2(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  ComplexMatrix X(SZ, SZ);
  ComplexMatrix Y(SZ, SZ);

  octave_value_list L;
  octave_value_list result;

  L.append(SZ);
  L.append(SZ);
  result = Frandn(L, L.length());
  X = result(0).complex_array_value ();  
  std::cout << "X\n"; 
  std::cout << X; 

  L.clear();
  L.append(X);
  result = Ffft2(L, L.length());
  Y = result(0).complex_array_value ();  
  std::cout << "Y\n"; 
  std::cout << Y; 

  return 0;
}

◆ 実行結果の例

Convex Hull (DLD_FUNCTIONS の convhulln)

convhull.a の中の Fconvhull() を使う。

◆ ソースコードの例

#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/ov-flt-complex.h>
#include<octave/ov-complex.h>
#include<octave/fCMatrix.h>
#include<octave/CMatrix.h>
#include<octave/defun-dld.h>

#define DIM 3
#define NPOINTS 40

extern octave_value_list Frandn(octave_value_list const&, int); 
extern octave_value_list Fconvhulln(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(NPOINTS, DIM);
  Matrix Y(NPOINTS, DIM);

  octave_value_list L;
  octave_value_list result;

  L.append(NPOINTS);
  L.append(DIM);
  result = Frandn(L, L.length());
  X = result(0).array_value ();  
  std::cout << "X\n"; 
  std::cout << X; 

  L.clear();
  L.append(X);
  result = Fconvhulln(L, L.length());
  Y = result(0).array_value ();  
  std::cout << "Y\n"; 
  std::cout << Y; 

  return 0;
}

◆ 実行結果の例

経過時間計測の例

ここでは、Octave の関数 Fconv2() を用いて， 2次元の畳み込み（コンボリューション）を行う．

【要点】

関数 conv2 は，Octave の conv2.a に実装されているので，「extern octave_value_list Fconv2(octave_value_list const&, int);」を含める．
「clock_t c = clock();」と「fprintf( stderr, "done, elapsed time = %f [sec]", ( (double)clock() - (double)c ) / CLOCKS_PER_SEC ); 」の2行は経過時間の測定のためです．

// liboctave を用いた 2次元畳み込み（コンボリューション）
#include<stdio.h>
#include<time.h>
#include<iostream>
#include<octave/config.h>
#include<octave/Matrix.h>
#include<octave/dMatrix.h>
#include<octave/defun-dld.h>

#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#include<gsl/gsl_statistics.h>


extern octave_value_list Fconv2(octave_value_list const&, int); 


// 正方行列のサイズを指定
#define DIM 3
#define DIM2 2

void matrix_rng_uniform( Matrix& M, const int N )
{
  // ここからは，GSL の機能を使って，乱数を配列 X, B に格納
  gsl_rng_env_setup();
  gsl_rng_type *T = (gsl_rng_type *)gsl_rng_default;
  /* 乱数発生器 */
  gsl_rng *r = gsl_rng_alloc(T);
  /* システムクロックを使って乱数の初期値を設定 */
  gsl_rng_set (r, time (NULL));

  int i, j;
  for(i = 0; i < N; i++) {
    for(j = 0; j < N; j++) {
      M(i, j) =  gsl_rng_uniform(r); 
    }
  }

  return;
}

int main( int argc, char **argv )
{

  Matrix X(DIM, DIM);
  Matrix B(DIM2, DIM2);
  Matrix Z(DIM, DIM);

  std::cout << X;
  std::cout << B;

  matrix_rng_uniform( X, DIM );
  matrix_rng_uniform( B, DIM2 );

  fprintf( stderr, "start, \n" );
  fprintf( stderr, "Z = Fconv2(X, B, \"full\"\n" );
  clock_t c = clock();

  octave_value_list L;
  L.append(X);
  L.append(B);
  L.append("full");
  octave_value_list result = Fconv2(L, L.length());
  Z = result(0).float_array_value (); 
  std::cout << Z; 

  fprintf( stderr, "done, elapsed time = %f [sec]\n", ( (double)clock() - (double)c ) / CLOCKS_PER_SEC );

  return 0;
}

ビルドの手順

Ubuntu で Intel Math Kernel Library(インテルマス・カーネル・ライブラリー) をインストール済みで，それを使いたい場合

# Linux, Intel MKL
g++ -O3 -I/usr/local/include/octave-3.6.2 -c -o oct.o oct.cc
g++ -O3 -o a.out oct.o -L/usr/local/lib/octave/3.6.2 -lgsl -loctave -lcruft -loctinterp -lreadline -L/opt/intel/lib/ia32 -L/opt/intel/mkl/lib/ia32 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread

Ubuntu で，ATLAS をパッケージを用いてインストール済みで，それを使いたい場合

g++ -O3 -I/usr/local/include/octave-3.6.2 -c -o oct.o oct.cc
g++ -O3 -o a.out oct.o -L/usr/local/lib/octave/3.6.2 -lgsl -loctave -lcruft -loctinterp -lreadline -L/usr/lib/atlas-base -L/usr/local/lib -llapack_atlas -llapack -lptf77blas -latlas -lmetis -lgfortran -lgomp -lpthread -ldl

Ubuntu で，ATLAS をソースコードからビルドしてインストール済みで，それを使いたい場合

g++ -O3 -I/usr/local/include/octave-3.6.2 -c -o oct.o oct.cc
g++ -O3 -o a.out oct.o -L/usr/local/lib/octave/3.6.2 -lgsl -loctave -lcruft -loctinterp -lreadline -L/usr/local/lib -L/usr/local/atlas/lib -ltmglib -llapack -lptcblas -lptf77blas -latlas -lmetis -lgfortran -lgomp -lpthread -ldl

Windows, Cygwin, ATLAS

-lstdc++ -lsupc++を付け忘れると，コンパイラのバージョンによってはエラーが出ることがあります．

g++ -I/usr/local/include/octave-3.6.2 -mthreads -c -o oct.o oct.cc
g++ --o a.out oct.o -mthreads -Xlinker --enable-auto-import -L/usr/local/lib/octave/3.6.2 -lgsl -loctave -lcruft -loctinterp -lreadline -ldl -L/usr/atlas/lib -lptf77blas -latlas -llapack_LINUX -ltmglib_LINUX -lpthread -lstdc++ -lsupc++

◆ 実行結果の例

OpenMP を使うように DLD-FUNCTIONS/conv2.cc を書き換える．2 次元の畳み込みを高速

【要点】

要するに，OpenMP 用に conv.cc を書き換えた（ソースコードは，このリンクからダウンロードできる）．コンパイルは通りました．正しい結果が出ているかのチェックは出来ていません．
g++ を使う．コンパイル時に「-fopenmp」を付ける．
4並列 (4スレッド) を前提にする　（これ以上スレッドを増やすと，書き換えがしんどい．主観ですが）．
conv2.cc の中のメソッド conv2 選び書き換える．

for ループを，等価な 4 つの for ループに分割．omp parallel, omp sections, omp section を使う．

#ifdef _OPENMP
#pragma omp parallel num_threads(4)
#endif
#ifdef _OPENMP
#pragma omp sections
#endif
{
#ifdef _OPENMP
#pragma omp section
#endif
  for (octave_idx_type oi1 = 0; oi1 < (outM / 4); oi1++)

...

#ifdef _OPENMP
#pragma omp section
#endif
  for (octave_idx_type oi2 = (outM / 4); oi2 < (outM / 2); oi2++)

...

#ifdef _OPENMP
#pragma omp section
#endif
  for (octave_idx_type oi3 = (outM / 2); oi3 < (outM / 1.33); oi3++)

...

#ifdef _OPENMP
#pragma omp section
#endif
  for (octave_idx_type oi4 = (outM / 1.33); oi4 < outM; oi4++)

for ループの中の局所変数 oi, oj sum, ai, bi, aj, bj, Ad, Bd は，oi1, oi2, oi3, oi4, oj1, ... のように名前を付け替える.（別スレッド間で同じ局所変数を共有しない）

やはり，並列化に，それなりの手間はかかる．やはり，箇所を厳選して書き換えるという方針が良い．
とは言え，OpenMP は便利．

Eclipse を使う場合

Eclipse のインストールが済んでいること．

「Octave の機能を呼び出す C++ プログラムの見本と実行手順」の Web ページの記述に従って，プロジェクトの作成と，プロジェクトのプロパティーの設定が済んでいること．

あらかじめ決めておく事項

プロジェクト名： Hello
プロジェクト・タイプ： 実行可能

Eclipse でのビルドと実行の手順

（前準備）「Octave の機能を呼び出す C++ プログラムの見本と実行手順」の Web ページの記述に従って，プロジェクトの作成と，プロジェクトのプロパティーの設定が済んでいること．
すべてビルド
CTRL + B（コントロールキーを押しながら「B」）または，「プロジェクト」メニューで「すべてビルド」
「コンソール」にコンパイルコマンド列が表示されるので，エラーが無いことを確認しておく．
実行
デフォルトでは，実行ファイルは， C:\workspace\Hello\Debug\Hello.exe に出来る．（「Hello」はプロジェクト名）
Eclipse のプロジェクト・エクスプローラで「Debug」を展開．Hello.exe を右クリックし，「実行」から「ローカル C/C++ アプリケーション」を選ぶ．数秒ほど待つと，コンソールに結果が出る．
※ Windows の Cygwin のコンソールを使うこともできる． C:\workspace\Hello\Debug\Hello.exe を実行．数秒ほど待つと，コンソールに結果が出てくる．
※ Cygwin で実行がうまくいかない場合， C:\cygwin\bin\cygwin1.dll を C:\Windows\System32 にコピーすることを忘れている可能性がある．

（参考）liboctave のヘッダファイルのうち関係部分の抜粋（Octave バージョン3.2.0）

ソースコードの引用に先立つ著作権表示

/*

Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007
              John W. Eaton
Copyright (C) 2008, 2009 Jaroslav Hajek

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, see
.

*/

Array.h

Array クラスは配列． MArray2 クラスの親クラスが Array2 クラスで，Array2 クラスの親クラスが Array

template <class T>
class
Array
{
  /* 省略 */
  void fill (const T& val);
  /* 省略 */
  octave_idx_type rows (void) const { return dim1 (); }
  octave_idx_type cols (void) const { return dim2 (); }
  /* 省略 */
  const dim_vector& dims (void) const { return dimensions; }
  /* 省略 */
  T& elem (octave_idx_type i, octave_idx_type j) { return elem (dim1()*j+i); }
  /* 省略 */
  T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); }
  /* 省略 */
  void print_info (std::ostream& os, const std::string& prefix) const;
};

MArray2.h

Marray2クラスは 2 次元配列．Matrix クラスの親クラス

template <class T>
class
MArray2 : public Array2<T>
{
  /* 省略 */
  MArray2 (octave_idx_type n, octave_idx_type m) : Array2<T> (n, m) { }
  MArray2 (octave_idx_type n, octave_idx_type m, const T& val) : Array2<T> (n, m, val) { }
  /* 省略 */
  MArray2<T>& insert (const Array2<T>& a, octave_idx_type r, octave_idx_type c)
  {
    Array2<T>::insert (a, r, c);
    return *this;
  }
};

dMatrix.h

Matrix クラスはdouble 型の 2 次元配列．

class
OCTAVE_API
Matrix : public MArray2<double>
{
  /* 省略 */
  Matrix (octave_idx_type r, octave_idx_type c) : MArray2<double> (r, c) { }
  Matrix (octave_idx_type r, octave_idx_type c, double val) : MArray2<double> (r, c, val) { }
  /* 省略 */
  Matrix& insert (const Matrix& a, octave_idx_type r, octave_idx_type c);
  Matrix& insert (const RowVector& a, octave_idx_type r, octave_idx_type c);
  Matrix& insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c);
  Matrix& insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c);
  /* 省略 */
  friend OCTAVE_API Matrix real (const ComplexMatrix& a);
  friend OCTAVE_API Matrix imag (const ComplexMatrix& a);
  /* 省略 */
  Matrix transpose (void) const { return MArray2<double>::transpose (); }
  /* 省略 */
  RowVector row (octave_idx_type i) const;
  ColumnVector column (octave_idx_type i) const;
  /* 省略 */
  Matrix inverse (void) const;
  /* 省略 */
  ComplexMatrix fourier2d (void) const;
  ComplexMatrix ifourier2d (void) const;
  /* 省略 */
  // Generic interface to solver with probing of type
  Matrix solve (const Matrix& b) const;
  ComplexMatrix solve (const ComplexMatrix& b) const;
  ColumnVector solve (const ColumnVector& b) const;

EIG.h

class
OCTAVE_API
EIG
{
  /* 省略 */
  EIG (const Matrix& a, bool calc_eigenvectors = true)
    { init (a, calc_eigenvectors); }
  /* 省略 */
  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
    { init (a, calc_eigenvectors); }
  /* 省略 */
  ComplexColumnVector eigenvalues (void) const { return lambda; }
  ComplexMatrix eigenvectors (void) const { return v; }
  /* 省略 */
};

dbleSVD.h

class
OCTAVE_API
SVD
{
  /* 省略 */
  enum type
    {
      std,
      economy,
      sigma_only
    };
  /* 省略 */
  SVD (const Matrix& a, type svd_type = SVD::std) { init (a, svd_type); }
  /* 省略 */
  DiagMatrix singular_values (void) const { return sigma; }
  Matrix left_singular_matrix (void) const;
  Matrix right_singular_matrix (void) const;
  /* 省略 */
};

本サイトは金子邦彦研究室のWebページである．

資料等の公開では，原則，「クリエイティブコモンズ BY NC SA」として公開するようにしている． PDFファイル，パワーポイントファイルなどには，「クリエイティブコモンズ BY NC SA」を明記するとともに，ロゴを記載するようにしている（作業が間に合っていない分もあるのでご容赦ください）．

公開している資料をご利用になる場合の，再配布の条件，剽窃の防止などについて，別ページ »で説明再配布や資料改変の際には，そのページをご確認ください．

サイトマップは，サイトマップのページをご覧下さい．本サイト内の検索は，サイト内検索のページをご利用下さい．

問い合わせ先：金子邦彦（かねこくにひこ） [image]