GSL のベクトル,行列に関する機能について,一部分を紹介する.
GNU Scientific Library -- Reference Manual の 20 章 「Statistics」
GNU Scientific Library -- Reference Manual の 17 章 「Random Number Generation」
GNU Scientific Library -- Reference Manual の 19 章 「Random Number Distributions」
GNU Scientific Library -- Reference Manual の 8 章 「Vectors and Matrices」
GNU Scientific Library -- Reference Manual の 12 章 「BLAS Support」
GNU Scientific Library -- Reference Manual の 14 章 「Eigensystems」
Windows での GSL のインストール(ソースコードを使用)(MSYS2,configure,make を利用): 別ページ »で説明している.
【要点】
#include<gsl/gsl_statistics_double.h>
飛び幅を 1 より大きい値に設定すると、double の配列の中の値を、等間隔で飛び飛びに使う。
スタートメニューで,「MSYS2」の下の「MSYS2 MINGW64」を選ぶ.
#include<stdio.h> #include<gsl/gsl_statistics.h> #define LEN 6 int main(int argc, char** argv) { double data[LEN] = {10.5, 18.2, 10.3, 15.4, 16.2, 18.3}; double max = gsl_stats_max( data, 1, LEN ); double min = gsl_stats_min( data, 1, LEN ); double mean = gsl_stats_mean( data, 1, LEN ); double sd = gsl_stats_sd( data, 1, LEN ); printf( "max: \t%f \n", max ); printf( "min: \t%f \n", min ); printf( "mean: \t%f \n", mean); printf( "sd: \t%f \n", sd ); return 0; }
gcc -I"C:\gsl\include" -c -o a.o a.c gcc -L"C:\gsl\lib" -o a.exe a.o -lgsl -lgslcblas ./a.exe
データ長 N で,飛び幅 1 のデータセット a の最大値.
gsl_stats_max(a, 1, N)
データ長 N で,飛び幅 1 のデータセット a の最小値.
gsl_stats_min(a, 1, N)
データ長 N で,飛び幅 1 のデータセット a を照準に整列.
gsl_sort(a, 1, N);
データ長 N で,飛び幅 1 のデータセット a を照準に整列. 置換を p に入れる.p は,置換を保持するのに十分な大きさ(つまり,サイズ N)の size_t 型の配列.
元の配列 a は変化しない
gsl_sort(p, a, 1, N);
上位3個を除き,0 で置き換え.
gsl_sort(p, a, 1, N); for( i = 0; i < N - 3; i++) { a[p[i]] = 0.0; }
【関連する外部ページ】
【要点】
#include<gsl/gsl_rng.h> #include<gsl/gsl_randist.h> #include<gsl/gsl_statistics.h>
スタートメニューで,「MSYS2」の下の「MSYS2 MINGW64」を選ぶ.
#include<stdio.h> // この3行は必須 #include<gsl/gsl_rng.h> #include<gsl/gsl_randist.h> #include<gsl/gsl_statistics.h> #define LEN 100 int main(int argc, char** argv) { double a[LEN]; double b[LEN]; 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; for(i=0; i<LEN; i++) { a[i] = gsl_rng_uniform(r); } for(i=0; i<LEN; i++) { b[i] = gsl_ran_gaussian(r, 1.0); } for(i=0; i<LEN; i++) { printf( "%d, \t%7.5f, \t%7.5f \n", i, a[i], b[i] ); } printf( "--, \t-------, \t%-------\n") ; printf( "mean: \t%7.5f, \t%7.5f \n", gsl_stats_mean(a, 1, LEN), gsl_stats_mean(b, 1, LEN) ); printf( "SD: \t%7.5f, \t%7.5f \n", gsl_stats_sd(a, 1, LEN), gsl_stats_sd(b, 1, LEN) ); gsl_rng_free(r); }
gcc -I"C:\gsl\include" -c -o a.o a.c gcc -L"C:\gsl\lib" -o a.exe a.o -lgsl -lgslcblas ./a.exe
gsl_rng_type *T = (gsl_rng_type *)gsl_rng_default; /* 乱数発生器 */ gsl_rng *r = gsl_rng_alloc(T);
上のように書いているので、環境変数 GSL_RNG_TYPE が設定されている場合には,それが乱数発生器として使われる. 環境変数 GSL_RNG_TYPE が未定義の場合,乱数発生器として mt19937 (メルセンヌ・ツイスター)が指定されたものとして振舞う.
解放は、gsl_rng_free()関数を使う。
gsl_rng_set (r, time (NULL));
time() 関数を使って設定している。
乱数発生器 r を使って,範囲 [0.0, 1.0) の一様分布の乱数を発生.
gsl_rng_uniform(r)
乱数発生器 r を使って,σ=s のガウス分布の乱数を発生.
gsl_ran_gaussian(r, s)
【要点】
ベクトルは,同じデータ型のデータの並び. 変数 v に,double 型の要素 2, 3, 4, 1 のベクトル(要素がこの順で並ぶ)を格納したいときは, 次のように書く.
gsl_vector *v = gsl_vector_alloc(4); // gsl_vector_alloc() は,double 型の要素を持つ GSL ベクトルの確保.中身は初期化されない. gsl_vector_set(v, 0, 2.0); gsl_vector_set(v, 1, 3.0); gsl_vector_set(v, 2, 4.0); gsl_vector_set(v, 3, 1.0); gsl_vector_free(v); // 領域の解放
下の1行は必須
#include<gsl/gsl_vector.h>
gsl_vector_double *v = gsl_vector_alloc(4); // gsl_vector_alloc() は,double 型の要素を持つ GSL ベクトルの確保.中身は初期化されない. gsl_vector_free(v); // 領域の解放
ベクトル v の第 i 要素を 1.0 にセット.i の値は,0 から,<ベクトルの長さ−1>.
gsl_vector_set( x, i, 1.0 );
ベクトルv の全要素を x にセット.
gsl_vector_set_all( v, x );
ベクトル src を,ベクトル dest にコピー
gsl_vector_memcpy( src, dest );
ベクトルv の第 i 要素を取り出す.i の値は,0 から,<ベクトルの長さ−1>.
double d = gsl_vector_get( v, i );
ベクトル b の要素の値を,ベクトル a の要素に加える (ai <- ai + bi)
gsl_vector_add( a, b )
ベクトル b の要素の値を,ベクトル a の要素から減ずる (ai <- ai - bi)
gsl_vector_sub( a, b )
ベクトル b の要素の値を,ベクトル a の要素に乗ずる (ai <- ai * bi)
gsl_vector_mul( a, b )
ベクトル b の要素の値を,ベクトル a の要素に除ずる (ai <- ai / bi)
gsl_vector_div( a, b )
ベクトル v の要素の値を x 倍する (vi <- x * vi)
gsl_vector_scale( v, x )
ベクトル v の最大値を得る
double d = gsl_vector_max( v );
ベクトル v の最大値を得る
double d = gsl_vector_min( v );
標準出力に,ベクトル v をバイナリ形式で書き込む
if ( gsl_vector_fwrite( stdout, v ) != 0 ) { fprintf( "error \n" ); }
標準入力から,ベクトル v をバイナリ形式で読み出す.
if ( gsl_vector_fread( stdin, v ) != 0 ) { fprintf( "error \n" ); }
標準出力に,ベクトル v をテキスト形式で書き込む
if ( gsl_vector_fprintf( stdout, v, "%f, " ) != 0 ) { fprintf( "error \n" ); }
標準入力から,ベクトル v をテキスト形式で読み出す.
if ( gsl_vector_fscanf( stdin, v, "%f, " ) != 0 ) { fprintf( "error \n" ); }
【関連する外部ページ】
【要点】
変数 M に,要素 1, 2, 3, 4 が入った2行2列の行列,つまり, 1 3 2 4 を格納したいときは,次のように書く. 「gsl_matrix_alloc(2, 2)」の最初の 2 が行数.次の 2 が列数.
gsl_matrix *M = gsl_matrix_alloc(2, 2); // gsl_matrix_alloc() は,double 型の要素の確保.中身は初期化されない. gsl_matrix_set(M, 0, 0, 1.0); gsl_matrix_set(M, 1, 0, 2.0); gsl_matrix_set(M, 0, 1, 3.0); gsl_matrix_set(M, 1, 1, 4.0); gsl_matrix_free(M); // 領域の解放
または
double d = { 1.0, 3.0, 2.0, 4.0 }; gsl_matrix_view MV = gsl_matrix_view_array( a, 2, 2 ); gsl_matrix v = &MV.matrix
下の1行は必須
#include<gsl/gsl_matrix.h>
gsl_matrix *M = gsl_matrix_alloc(2, 2); // gsl_matrix_alloc() は,double 型の要素の確保.中身は初期化されない. gsl_matrix_free(M); // 領域の解放
行列 M の第 i 行, 第 j 列の要素を 1.0 にセット.i, j の値は,0 から,行列の縦,横−1.
gsl_matrix_set( M, i, j, 1.0 );
行列 の全要素を x にセット.
gsl_matrix_set_all( M, x );
つまり,対角要素を 1 に,非対角要素を 0 にセットする.行列 M が正方行列でなくても動く.
gsl_matrix_set_identity( M );
行列 src を,行列 dest にコピー
gsl_matrix_memcpy( dest, src );
ベクトル像は,ベクトルと同様に扱うことができる.
行列 M の第 j 行を取り出す.これは参照を使ってベクトル像を生成している.
※ gsl_matrix v = &v1.matrix のように書いて,gsl_martix_double へのポインタが得られる
gsl_vector_view v1; v1 = gsl_matrix_row(T, i)
行列 M の第 j 列を取り出す.これは参照を使ってベクトル像を生成している.
※ gsl_matrix v = &v1.matrix のように書いて,gsl_martix_double へのポインタが得られる
gsl_vector_view v1; v1 = gsl_matrix_column(T, i)
行列 M の i 行 j 列の要素を取り出す.i の値は,0 から,<行列の行数−1>.j の値は,0 から,<行列の列−1>.
double d = gsl_matrix_get( v, i, j );
標準出力に,行列 M をバイナリ形式で書き込む
if ( gsl_matrix_fwrite( stdout, v ) != 0 ) { fprintf( "error \n" ); }
標準入力から,行列 M をバイナリ形式で読み出す.
if ( gsl_matrix_fread( stdin, v ) != 0 ) { fprintf( "error \n" ); }
標準出力に,行列 M をテキスト形式で書き込む
if ( gsl_matrix_fprintf( stdout, v, "%f, " ) != 0 ) { fprintf( "error \n" ); }
標準入力から,行列 M をテキスト形式で読み出す.
if ( gsl_matrix_fscanf( stdin, v, "%f, " ) != 0 ) { fprintf( "error \n" ); }
【関連する外部ページ】
BLAS(Basic Linear Algebra Subprograms)とは,行列演算,ベクトル演算の機能(下記)をもったプログラム群. GSL には,BLAS の機能を呼び出す機能があります.
ベクトルとベクトルの演算
DOT : 内積, AXPY : AXPY 演算 ( y <- ax + y の形など), NORM : ノルム など
行列とベクトルと計算
行列とベクトルの積 ( y <- Ax ), 行列の rank-1 更新 ( A <- A + xy' ) など
行列同士の演算
行列と行列の積 ( Z <- XY ) など
GSL ベクトル v1 と v2 の内積(ドット積). 「ddot」の「d」は,GSL ベクトル v1 と v2 が double の精度という意味.
#include<gsl_blas.h> ... double result; gsl_blas_ddot(v1, v2, &result)
GSL ベクトルv のノルム(ユークリッド距離によるノルム; Euclidean norm). 「dnrm2」の「d」は,GSL ベクトル v が double の精度という意味.
#include<gsl_blas.h> ... double n = gsl_blas_dnrm2(v);
GSL ベクトル v の要素の総和. 「dasum」の「d」は,GSL ベクトル v が double の精度という意味.
#include<gsl_blas.h> ... double n = gsl_blas_dasum(v);
GSL ベクトル v の要素を alpha 倍して,v 二格納. 「dscal」の「d」は,GSL ベクトル v が double の精度という意味.
#include<gsl_blas.h> ... double n = gsl_blas_dscal( alpha, v );
y = M x + y
#include<gsl_blas.h> ... /* y = M x + y */ gsl_blas_dgemv( CblasNoTrans, 1.0, M, x, 1.0, y );
C = AB + C
#include<gsl_blas.h> ... /* C = A B */ gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C );
【関連する外部ページ】:
【要点】
#include<gsl/gsl_eigen.h>
D × D の実数対称行列の固有値を計算するための作業領域の確保
gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc( D );
GSL 行列 (gsl_matrix) A の固有値を計算する. A は実数対称行列であること. 計算された固有値は,GSL ベクトル (gsl_vector) v に格納される.
実行の結果 A は破壊される.
gsl_eigen_symm( A, v, w );
D × D の実数対称行列の固有値と固有ベクトルを計算するための作業領域の確保
gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc( D );
GSL 行列 (gsl_matrix) A の固有値と固有ベクトルを計算する. A は実数対称行列であること. 計算された固有値は,GSL ベクトル (gsl_vector) v に格納される. 計算された固有ベクトルは,GSL 行列 (gsl_matrix) c に列として格納される. v の要素と,c の列は対応する(例えば,vの第一要素は,c の第一列になる). 固有ベクトルは,お互いに直交し,長さが 1 になるように正規化されている.
実行の結果 A は破壊される.
gsl_eigen_symmv( A, v, c, w );
GSL ベクトル (gsl_vector) v に格納されている固有値と, GSL 行列 (gsl_matrix) c 格納されている固有ベクトルを整列する.
gsl_eigen_symmv_sort( v, c, GSL_EIGEN_SORT_ABS_DEC );
作業領域の解放
gsl_eigne_symmv_free( w );
D × D の実数非対称行列の固有値を計算するための作業領域の確保
gsl_eigen_nonsymm_workspace *w = gsl_eigen_nonsymm_alloc( D );
GSL 行列 (gsl_matrix) A の固有値を計算する. A は実数非対称行列. 計算された固有値は,GSL 複素数ベクトル (gsl_complex_vector) v に格納される.
実行の結果 A は破壊される.
gsl_eigen_nonsymm( A, v, w );
D × D の実数非対称行列の固有値と固有ベクトルを計算するための作業領域の確保
gsl_eigen_nonsymmv_workspace *w = gsl_eigen_nonsymmv_alloc( D );
GSL 行列 (gsl_matrix) A の固有値と固有ベクトルを計算する. A は実数非対称行列であること. 計算された固有値は,GSL 複素数ベクトル (gsl_complex_vector) v に格納される. 計算された固有ベクトルは,GSL 複素数行列 (gsl_complex_matrix) c に列として格納される. v の要素と,c の列は対応する(例えば,vの第一要素は,c の第一列になる). 固有ベクトルは,お互いに直交し,長さが 1 になるように正規化されている.
実行の結果 A は破壊される.
gsl_eigen_nonsymmv( A, v, c, w );
GSL ベクトル (gsl_vector) v に格納されている固有値と, GSL 行列 (gsl_matrix) c 格納されている固有ベクトルを整列する.
gsl_eigen_nonsymmv_sort( v, c, GSL_EIGEN_SORT_ABS_DEC );
作業領域の解放
gsl_eigne_nonsymmv_free( w );
【関連する外部ページ】