数値計算マニュアル：線型代数 (1)

更新履歴と目次

更新履歴

2008-10-13

• LD? (Leading Dimension of ?)に関する注意事項を追加しました．
• サンプルコードを更新しました．

2008-03-05

サンプルコード中のMakefileを修正し，FreeBSD 6.3Rに対応しました．

2005-10-21

サンプルコードを修正しました．

2005-08-23

サンプルコード中のMakefileを修正しました．

2005-07-23

ルーチンcblas_?syr2，cblas_?trmm，cblas_?syr2kの解説を追加しました． これに合わせてサンプルコードを更新しました． 実行列の一般的な演算はほとんど網羅されました． 行列の型に配慮することで，複素系を除くほぼ全てのルーチンを利用可能なはずです．

2005-07-20

ルーチンcblas_?trsmの解説を追加しました． これに合わせてサンプルコードを更新しました．

2005-06-28

ルーチンcblas_?symv，cblas_?syr，cblas_?trmv，cblas_?trsv，cblas_?syrkの解説を追加しました． これに合わせてサンプルコードを更新しました．

目次

• はじめに
• 利用の前に
• Level 1 BLAS
• Level 2 BLAS
• Level 3 BLAS
• ライブラリのリンクとサンプルコード
• おわりに

はじめに

BLASとは

BLASはBasic Linear Algebra Subprogramsの略で線型代数の基本的な演算を行なうためのライブラリ群です． この中には行列の積や和，行列とベクタの積などの行列計算が含まれます． Fortranで開発されており，数値計算用ライブラリとして精度，安定性に定評があります． ここではBLASをC言語から利用するためのインターフェース，CBLASについて説明します． このマニュアルで扱うCBLASはATLASによって提供されるC言語インターフェースを想定していますが， intelのMKL (Math Kernal Library)やGSL (GNU Scientific Library) のCBLASもほぼ同じ方法で使うことができると思われます． これはCBLASにはインターフェース規約というものが存在し， 関数の仕様が決められているためです．

BLASとATLAS

ATLASはAutomatically Tuned Linear Algebra Softwareの略でコンパイル時に自動的にCPUを始めとする， 計算機の実行条件に合わせて最適化を行なったBLASの一種です． SSE等のCPU特殊命令についても使用／不使用が自動的に設定されます． このため，一般にBLASよりも高速であり，通常の場合はこちらを使うことが推奨されます．

上述のようにATLASはBLASの上位互換と考えることができます． 加えて，ATLASにはLAPACKの一部のルーチンも同時に最適化を行ないます． これらの最適化ルーチンについてはライブラリのリンクのところで詳述します．

インストール

本研究室の標準インストールに従っていれば，BLAS，ATLASはインストールされているはずです， 何らかの理由でインストールされていない場合は管理者に相談してください． 一般の方はLAPACK (FAQ)を参照してください． ここで，Visual Studioと心中する覚悟のある方は， 番外編 (1)を参照してください．

利用の前に

利用の前に注意しなければならないことについて解説します．

Fortranの規則

行形式と列形式

二次元配列をメモリ上に格納する際のルールがFortranとC言語とでは異なります． 以下の行列(二次元配列)を例にとり説明します． Fortranは列形式でこれを格納します． つまりメモリアドレスの小さい順に， a，d，g，b，e，h，c，f，iとなります． C言語は行形式， a，b，c，d，e，f，g，h，iとなります． BLASのC言語インターフェースであるCBLASルーチンは引数によって行形式も列形式も扱えるようになっていますが， LAPACKとの連携を考えると列形式を利用すべきでしょう． 行列の形式の指定は列形式がCblasColMajor，行形式がCblasRowMajorです．

LD?に関する注意事項

本マニュアルはlda，ldb等，LD?なる変数を「行列の行数」として言及していますが， これは本来は「Leading Dimension of ?」，「配列?の第一次元」であり， メモリ空間上で連続的に配置される要素の要素数を意味します． 二次元配列の格納形式に列形式(CblasColMajor)を採用した場合， メモリ空間上に連続的に配置されるのは列要素ですから，その要素数たるLD?は「行列の行数」となります． しかしながら，C言語の標準の格納形式である行形式(CblasRowMajor)を採用した場合， メモリ空間上に連続的に配置されるのは行要素ですから， LD?は「行列の列数」となります． したがって，格納形式を変更する場合はLD?を適宜読み替えてください．

また，Fortranの二次元配列の要素は必ず連続領域でなければなりませんが， C言語は二次元配列は一次元配列へのポインタへの配列を意味しますから， 自動変数として確保した場合を除けば，要素の実体の領域が連続であるとは限りません． したがって，二次元配列A[i][j]などとして作った行列を， ポインタとしてCBLASの関数に与える場合は， 格納形式だけでなく必ず実体が連続的に配置されていることもまた確認してください．

引数の形式

Fortranはサブルーチンのコールの際に値渡しを行なうことができません． 引数は全てポインタとなります． しかしながら，CBLASルーチンはこのことが考慮されており，値で引数を渡せるようになっています．

変数の型と精度

変数の型と精度は以下の表のようになります．複素数は正味ポインタとなる，と考えてください． 各種ルーチンの?の部分に入ります．

このほかに，関数の場合は引数と戻り値の型が異なる場合があります． ds等と表記する場合がそれです，これは，単精度実数を引数として倍精度実数を返すことを意味します．

行列の型

CBLASは特殊な型の行列に対してより高速な演算を行なうルーチンを提供します． 以下の表にそうした行列の型を挙げます．

行列の型とその略号

名称	略号	例	説明
一般行列(General Matrix)	ge		要素の配置に制約のない一般的な行列です．
一般帯行列(General Band Matrix)	gb		非零の要素が主対角線上とそれに平行な帯の中に存在する行列です． 応用上良く現れます．左の例は3重対角行列です．
対称行列(Symmetric Matrix)	sy		転置行列がもとの行列と等しくなる行列です．主に実行列のことを指します． やはり，応用上良く現れます．
対称帯行列(Symmetric Band Matrix)	sb		転置行列がもとの行列と等しくなる帯行列です．
エルミート行列(Hermitian Matrix)	he		転置共役行列がもとの行列と等しくなる行列です．主に複素行列ですが対角要素は実数です． 応用上良く現れます．
エルミート帯行列(Hermitian Band Matrix)	hb		転置共役行列がもとの行列と等しくなる帯行列です．
三角行列(Triangular Matrix)	tr		非零の要素が主対角線上の上のみ(下のみ)に存在する行列です． 左の例は上三角行列です．
三角帯行列(Triangular Band Matrix)	tb		非零の要素が主対角線上とその上の(下の)平行な帯の中に存在する行列です． 左の例は上三角帯行列です．

行列に対する事前の演算

CBLASが提供するルーチンの内，幾つかのものは行列に対する事前の演算を指定する必要があります． ここで用いる定数はcblas.h内で定義されています．

転置演算

与えられた行列の転置を演算に供します．引数の型はCBLAS_TRANSPOSE型です． 転置なし(CblasNoTrans)，転置(CblasTrans)，共役転置(CblasConjTrans)が存在します．

三角行列の演算

与えられた行列の上(下)三角部分のみを演算に供します．引数の型はCBLAS_UPLO型です． 上三角部分(CblasUpper)，下三角部分(CblasLower)が存在します． これにより，三角行列関係のルーチンは与えられた行列が上三角行列なのか下三角行列なのかを判断します． また，対称行列関係のルーチンも指定された部分のみを演算に供することで演算を省略／高速化します．

三角行列の演算(2)

三角行列関係のルーチンにおいて，与えられた行列を単位三角行列とみなして演算に供します． 一般三角行列(CblasNonUnit)，単位三角行列(CblasUnit)が存在します．

演算の向き

cblas_?trmm()とcblas_?trsm()において， 与えられた行列を演算に供する際に左から掛けるのか右から掛けるのかを指定します． 引数の型はCBLAS_SIDE型です． 左から(CblasLeft)，右から(CblasRight)が存在します．

ルーチンの分類

BLASには演算の種類によって三種類の分類が存在します． 大まかに次のように分類されています．

Level 1 BLAS

ベクタ─ベクタ演算ルーチンを提供します．

Level 2 BLAS

ベクタ─行列演算ルーチンを提供します．

Level 3 BLAS

行列─行列演算ルーチンを提供します．

Level 1 BLAS

既に述べた通り，Level 1 BLASはベクタ─ベクタ演算ルーチンを提供します． BLASは高速化のために演算に2オペランド方式を採用しています． このため，演算に供されるベクタの内，必ず一方の値が失われることを覚えておいてください．

ルーチンの説明は提供される精度を挙げたあと， 最も良く利用すると思われる精度のルーチンについて引数と戻り値を解説します．

Level 1 BLASの全ルーチン

まずLevel 1 BLASの全ルーチンを表にします．

void cblas_dswap(const int N,
                 const double *X,
                 const int incX,
                 double *Y,
                 const int incY);

void cblas_dscal(const int N,
                 const double alpha,
                 double *X,
                 const int incX);

void cblas_dcopy(const int N,
                 const double *X,
                 const int incX,
                 double *Y,
                 const int incY);

void cblas_daxpy(const int N,
                 const double alpha,
                 const double *X,
                 const int incX,
                 double *Y,
                 const int incY);

double cblas_ddot(const int N,
                  const double *X,
                  const int incX,
                  const double *Y,
                  const int incY);

void cblas_zdotu_sub(const int N,
                     const void *X,
                     const int incX,
                     const void *Y,
                     const int incY,
                     void *dotu);

void cblas_zdotc_sub(const int N,
                     const void *X,
                     const int incX,
                     const void *Y,
                     const int incY,
                     void *dotu);

double cblas_dnrm2(const int N,
                   const double *X,
                   const int incX);

double cblas_dasum(const int N,
                   const double *X,
                   const int incX);

CBLAS_INDEX cblas_idamax(const int N,
                         const double *X,
                         const int incX);

Level 2 BLAS

Level 2 BLASはベクタ─行列演算ルーチンを提供します． Level 1 BLAS同様に左辺値の計算対象が右辺にも存在します． このため，演算に供されるベクタまたは行列の内，必ず一つの値が失われることを覚えておいてください．

本節で初めて行列が登場します．CBLASは行形式，列形式双方の行列を扱うことができますが， LAPACKとの併用を考えた場合，列形式での利用をお勧めします． また，演算に供する行列が既述の特殊な型に該当する場合，専用のルーチンを利用した方が 速度／精度の面で有利になります．

特殊な型の行列，対称(エルミート)行列(sy，he)と三角行列(tr) はCBLAS_UPLOで指定した三角部分以外は不定値(正確には更新されない)となりますので特に注意が必要です． 必ずサンプルコードに目を通してください．

ルーチンの説明は提供される精度を挙げたあと， 最も良く利用すると思われる精度のルーチンについて引数と戻り値を解説します．

Level 2 BLASの全ルーチン

Level 2 BLASの全ルーチンを表にします．

ルーチン名を注意深く見ると分かりますが， 同じ処理内容のルーチンが行列の型に合わせて複数存在します． 例えばcblas_?gemvは General Matrix に対する Matrix-vector productの意味です． 同じ処理内容のルーチンが帯行列，対称行列，三角行列向けに存在することが分かります．

void cblas_dgemv(const enum CBLAS_ORDER Order,
                 const enum CBLAS_TRANSPOSE TransA,
                 const int M,
                 const int N,
                 const double alpha,
                 const double *A,
                 const int lda,
                 const double *X,
                 const int incX,
                 const double beta,
                 double *Y,
                 const int incY);

void cblas_dger(const enum CBLAS_ORDER Order,
                const int M,
                const int N,
                const double alpha,
                const double *X,
                const int incX,
                const double *Y,
                const int incY,
                double *A,
                const int lda);

void cblas_dsymv(const enum CBLAS_ORDER Order,
                 const enum CBLAS_UPLO Uplo,
                 const int N,
                 const double alpha,
                 const double *A,
                 const int lda,
                 const double *X,
                 const int incX,
                 const double beta,
                 double *Y,
                 const int incY);

void cblas_dsyr(const enum CBLAS_ORDER Order,
                const enum CBLAS_UPLO Uplo,
                const int N,
                const double alpha,
                const double *X,
                const int incX,
                double *A,
                const int lda);

void cblas_dsyr2(const enum CBLAS_ORDER Order,
                 const enum CBLAS_UPLO Uplo,
                 const int N,
                 const double alpha,
                 const double *X,
                 const int incX,
                 const double *Y,
                 const int incY,
                 double *A,
                 const int lda);

void cblas_dtrmv(const enum CBLAS_ORDER Order,
                 const enum CBLAS_UPLO Uplo,
                 const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_DIAG Diag,
                 const int N,
                 const double *A,
                 const int lda,
                 double *X,
                 const int incX);

void cblas_dtrsv(const enum CBLAS_ORDER Order,
                 const enum CBLAS_UPLO Uplo,
                 const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_DIAG Diag,
                 const int N,
                 const double *A,
                 const int lda,
                 double *X,
                 const int incX);

Level 3 BLAS

Level 3 BLASは行列─行列演算ルーチンを提供します． 左辺値の計算対象が右辺にも存在するのはLevel 2 BLASと同様です． 演算に供されるベクタや行列の内，必ず一つの値が失われることを覚えておいてください．

ルーチンの説明は提供される精度を挙げたあと， 最も良く利用すると思われる精度のルーチンについて引数と戻り値を解説します．

Level 3 BLASの全ルーチン

Level 3 BLASの全ルーチンを表にします．

void cblas_dgemm(const enum CBLAS_ORDER Order,
                 const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB,
                 const int M,
                 const int N,
                 const int K,
                 const double alpha,
                 const double *A,
                 const int lda,
                 const double *B,
                 const int ldb,
                 const double beta,
                 double *C,
                 const int ldc);

void cblas_dsyrk(const enum CBLAS_ORDER Order,
                 const enum CBLAS_UPLO Uplo,
                 const enum CBLAS_TRANSPOSE Trans,
                 const int N,
                 const int K,
                 const double alpha,
                 const double *A,
                 const int lda,
                 const double beta,
                 double *C,
                 const int ldc);

void cblas_dsyr2k(const enum CBLAS_ORDER Order,
                  const enum CBLAS_UPLO Uplo,
                  const enum CBLAS_TRANSPOSE Trans,
                  const int N,
                  const int K,
                  const double alpha,
                  const double *A,
                  const int lda,
                  const double *B,
                  const int ldb,
                  const double beta,
                  double *C,
                  const int ldc);

void cblas_dtrmm(const enum CBLAS_ORDER Order,
                 const enum CBLAS_SIDE Side,
                 const enum CBLAS_UPLO Uplo,
                 const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_DIAG Diag,
                 const int M,
                 const int N,
                 const double alpha,
                 const double *A,
                 const int lda,
                 double *B,
                 const int ldb);

void cblas_dtrsm(const enum CBLAS_ORDER Order,
                 const enum CBLAS_SIDE Side,
                 const enum CBLAS_UPLO Uplo,
                 const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_DIAG Diag,
                 const int M,
                 const int N,
                 const double alpha,
                 const double *A,
                 const int lda,
                 double *B,
                 const int ldb);

ライブラリのリンク法とサンプルコード

ライブラリのリンク

ライブラリの詳細を以下に示します．Lapack関係のライブラリもついでに挙げておきます．

以上をまとめると，ATLASを用いる場合はライブラリのリンクを以下のように行なえば良いということになります． FreeBSD(上)，Debian GNU/Linux(下)です．

-lalapack -lf77blas -lcblas -latlas -lf2c -lm

-llapack_atlas -llapack -lf77blas -lcblas -latlas -lf2c -lm

しかしながら，ATLASはコンパイル時にgccに対して-ffast-mathオプションを指定することがあるようです． ATLASを用いたところ，精度が悪化した，等のことがあった場合は汎用のBLAS，LAPACKを利用すべきかもしれません． 実際，FreeBSDではリンクするライブラリによって計算結果が異なる場合があります．

サンプルコード

サンプルコードを置いておきます．適当なところで展開するとMakefileと以下のソースコードが展開されます．

• general.c
• special.c
• disp_array.c

ここで，disp_array.cは行列の表示のためのモジュールです． general.cは一般行列，ベクタに関するサンプルコードです． special.cは対称行列，三角行列の演算に関するサンプルコードです． 詳しくはコード中のコメントを読んでください．

おわりに

本マニュアルは不完全です，改訂を歓迎します．

• 全ルーチンの演算内容を和訳する
• 複素行列についてマニュアルを補う

© Copyright 2006-2009 ej-sib All rights reserved.