数値計算第 4 回課題 (大石先生)

{G99P043-4 \atop 河邊昌彦}

出題日 : 2001/06/06
提出期限 : 2001/06/20
提出日 : 2001/06/20

1 問題

(1) 行列の積のプログラムを使って dgemm_を作れ。 CLAPACK をインストールしながら仕様を見よ。

(2) Ax = b を解くための LAPACK の関数

         dgesv_( &x, &inc, A, &n, piv, X, &n, &info )

を作れ (BLAS を高速化したら高速になる形で実装せよ)。

2 dgemm

2.1 仕様

?gemm はスカラー-行列-行列の積を、スカラー-行列の積の結果に足す関数である。具体的には、 op_{A}(A) を (m, k)-行列、 op_{B}(B) を (k, n)-行列、 C を (m, n)-行列とし、 α, β をスカラーとすると、

C := α · op_{A}(A) · op_{B}(B) + β · C

という計算をする。ただし、 op_{A}, op_{B} はそれぞれ、恒等写像か或いは転置である。そして、スカラーや行列の要素の型によって幾つかのバージョンがあり、倍精度浮動小数点数のものを扱う関数

が dgemm である。

C 言語における dgemm は次のような引数使用を持つ。

     int dgemm_( char *transa, char *transb, integer *m, integer *n,
                 integer *k, doublereal *alpha, doublereal *a, integer *lda,
                 doublereal *b, integer *ldb, doublereal *beta, doublereal *c,
                 integer *ldc );

ここで注意しなければならないのは、もともと標準 BLAS の関数は FORTRAN で書かれており、 f2c によって変換されている、ということである。従って、関数名の最後にはアンダースコアが付加され、引数は全てポインタ渡しになっている。

また、 FORTRAN の 2 次元配列の仕様により、行列の添字は 1 オリジンであり、同じ列の要素が連続したアドレスに格納されている。即ち、 FORTRAN における、行列の (i, j)-成分へのアクセス方法を C で書くと、

FORTRAN  C

      宣言  DOUBLE PRECISION A( LDA, * )  doublereal *a; integer lda;
      アクセス A(i,j)                        a[ i-1 + (j-1)*lda ]

となる。

なお、上で使われている doublereal, integer 等といった型は、 f2c.h で定義されている。ここで、 integer は long に typedef されているので、 int*型のポインタはそのままでは渡せないことに注意する。

個々のパラメータについては、次のようになっている。 *m, *n, *k は行列のサイズを表しており、前述のように a が*m 行*k 列、 b が*k 行*n 列、 c が*m 行*n 列である。 *alpha, *beta は行列に掛けるスカラーである。 a, b, c が行列であり、それぞれ*lda, *ldb, *ldc を使って、要素にアクセスする。

2.2 実装方法

FORTRAN と C で行列の格納方向が逆なので、 FORTRAN の (m, n)-行列をそのまま C の行列であると思ってアクセスすると、転置された (n, m)-行列に見えることになる。従って、 FORTRAN における行列 A_{F} , B_{F} , C_{F} に対して

α · op_{A}(A_{F}) · op_{B}(B_{F}) + β · C_{F}

という計算を行うと、 C の行列と見なした A_{C}, B_{C}, C_{C} に対しては、

α · op_{A}(A_{F}) · op_{B}(B_{F}) + β · C_{F} = α · op_{A}(^{t}A_{C}) · op_{B}(^{t}B_{C}) + β ·^{t}C_{C}
= α ·^{t}(op_{B}(B_{C}) · op_{A}(A_{C})) + β ·^{t}C_{C}
= ^{t}(α · op_{B}(B_{C}) · op_{A}(A_{C}) + β · C_{C})

という計算を行ったものと同じになる。よって、

C_{F} := α · op_{A}(A_{F}) · op_{B}(B_{F}) + β · C_{F}
⟺ C_{C} := α · op_{B}(B_{C}) · op_{A}(A_{C}) + β · C_{C}

ということである。

前回のレポートで作成した、行列の積を計算するプログラムは C の行列を扱うようになっている。そこで、今回作成する dgemm_の中では、

C_{C} := α · op_{B}(B_{C}) · op_{A}(A_{C}) + β · C_{C}

という計算をすることで、 FORTRAN の行列を扱うことにした。ただし、 op_{A}(A_{C}) は (k, m)-行列、 op_{B}(B_{C}) は (n, k)-行列、 C_{C} は (n, m)-行列になることに注意する必要がある。

2.3 プログラム

上のような方針で前回のプログラムを利用し、 dgemm_を実装した (my_dgemm) のが B.1 である。行列の大きさが前もって見積もれなくなったため、ブロック化は同じメモリ領域を使いまわすようにしたり、或いは、テンプレートを使って最適化しやすくしてみたりしたが、基本は前回のプログラムと同じである。また、スカラーや要素の型をテンプレートパラメータとしたので、すぐに sgemm_も実装することができる。

プログラムの検証は、既存の標準 BLAS 又は最適化 BLAS の dgemm_と結果を比べることで出来る。ここでは、 Intel の Math Kernel Library(以下 MKL) という Windows, PentiumIII 用の最適化 BLAS を使用した。

2.4 実行結果

実行結果を A.1 に記す。これに依れば、 my_dgemm の計算誤差には問題が無いといってよいであろう。また、倍精度浮動小数点の場合、 my_dgemm は最適化 BLAS に比べて 1.3〜1.5 倍程度の時間がかかるようである。一方、単精度浮動小数点の場合は、 PentiumIII にはストリーミング SIMD 拡張命令があるので、それを利用しているであろう最適化 BLAS の方が大幅に速い。

3 dgesv

3.1 仕様

  dgesv_は、

AX = B

という (複数の) 連立一次方程式を解くための関数である。ただし、 A は (n, n)-行列、 X, B は (n, n_{rhs})-行列である。これは、 X = [x_{1} x_{2} ⋯ x_{n_{rhs}}] , B = [b_{1} b_{2} ⋯ b_{n_{rhs}}] とすると、

Ax_{i} = b_{i}  (1 ≤ i ≤ n_{rhs})

という n_{rhs} 個の n 元連立一次方程式と見ることもできる。

dgesv_は次のように宣言されている。

     int dgesv_( integer *n, integer *nrhs, doublereal *a, integer *lda,
                integer *ipiv, doublereal *b, integer *ldb, integer *info )

入力である*n, *nrhs, a, *lda, b, *ldb は、上で示したものをそれぞれ表している。そして、 ipiv は配列であり、ピボット選択による行の交換が出力される。即ち、第 i 行は第 ipiv[i] 行と交換された、ということである。解 X は b の上に上書きされ、 a には元の a の LU 分解が上書きされる。 *info に結果の状態を表す値が入って返ってくるが、 *info が 0 ならば成功であり、正の値ならば i 列目のピボット選択に失敗したことを表す。また、負の値ならばパラメータの値が不正であったことを示している。

3.2 実装方法

上の

AX = B

において、 A が A = LU と LU 分解できたならば、

LUX = B

であるので、

\{ \begin{matrix}LY & = B \\ UX & = Y \end{matrix}

となる。片方の行列が三角行列である、行列の方程式を解くには、 BLAS の dtrsm_を使うことができる。従って、 A を LU 分解することができたならば、 dtrsm_によって X を求められる。

では、 BLAS を使って効率よく LU 分解するにはどのようにすればよいのだろうか。これは、行列の積と同様に、行列をブロック化して考えればよい。

A = [\begin{matrix}A_{11} & A_{12} & ⋯ & A_{1m} \\ A_{21} \\ ⋮ & & A′ \\ A_{m1} \end{matrix}]

とブロック化し、 LU 分解した結果も

L = [\begin{matrix}L_{1} & O & ⋯ & O \\ L_{2} \\ ⋮ & & L′ \\ L_{m} \end{matrix}]  U = [\begin{matrix}U_{1} & U_{2} & ⋯ & U_{m} \\ O \\ ⋮ & & U′ \\ O \end{matrix}]

とブロック化されたとする。ただし、 L_{1}, L′ は下三角行列、 U_{1}, U′ は上三角行列である。すると、

LU = [\begin{matrix}L_{1} & O & ⋯ & O \\ L_{2} \\ ⋮ & & L′ \\ L_{m} \end{matrix}] [\begin{matrix}U_{1} & U_{2} & ⋯ & U_{m} \\ O \\ ⋮ & & U′ \\ O \end{matrix}] = [\begin{matrix}A_{11} & A_{12} & ⋯ & A_{1m} \\ A_{21} \\ ⋮ & & A′ \\ A_{m1} \end{matrix}] = A

であるから、

[ L_{1} ; ⋮ ; L_{m}] U_{1} = [ A_{11} ; ⋮ ; A_{m1}]  (3.1)

L_{1} [U_{2} ⋯ U_{m}] = [A_{12} ⋯ A_{1m}]  (3.2)

[ L_{2} ; ⋮ ; L_{m}] [U_{2} ⋯ U_{m}] + L′U′ = A′  (3.3)

である。 (3.1) は、右辺の LU 分解が左辺となっている。また、 (3.2) は、 [U_{2} ⋯ U_{m}] を未知数とする方程式と見ることができる。そして、 (3.3) は、移項すれば

L′U′ = − [ L_{2} ; ⋮ ; L_{m}] [U_{2} ⋯ U_{m}] + A′ (3.4)

となり、右辺を LU 分解していることが分かる。

従って、まず、 (3.1) で普通に LU 分解をし、 L_{1}, ⋯, L_{m}, U_{1} を得る。次に、そこで得られた L_{1} を使い、 (3.2) において dtrsm_で U_{2}, ⋯, U_{m} を求める。そして、 dgemm_で(3.4) の右辺を計算する。後は、 dgemm_の結果を L′, U′ に LU 分解すればよいので、同様のことを繰り返す。

ここではピボット選択による行の交換は考慮しなかったが、考慮する場合でも同様である。よって、このようにして、 A の LU 分解を求めることができる。

3.3 プログラム

以上のような方針で my_dgesv を実装したプログラムが B.2 である。この中では、以下の BLAS の関数を使っている。

関数名機能
idamax_ベクトルの中で絶対値が最大の要素を探す。
dswap_二つのベクトルを交換する。
dscal_ベクトル × スカラーの計算。
dger_縦ベクトル × 横ベクトルの計算。
dtrsm_三角行列の方程式を解く。

                 dgemm_行列の積。

ここでは、 LU 分解のブロックのサイズは定数とした。しかし、実際の LAPACK における dgetrf_(LU 分解の関数) は、 dgemm_が高速に行われるように、 ilaenv_という関数で dgemm_のブロックのサイズを調べ、 LU 分解のブロックのサイズをそれに合わせるようにしている。

3.4 実行結果

標準 BLAS を使用した場合と、最適化 BLAS(MKL) を使用した場合の実行結果を A.2 に載せる。結果の検証は、 LAPACK(標準 BLAS 使用) の dgesv_と比較することで行った。行列 A, B の要素は毎回ランダムで決めたので、それに依存して実行時間が変わっている (ピボット選択にかかる時間が変わる) が、最適化 BLAS を使用すると標準 BLAS の 4〜8 倍の速さになった。

4 考察

今回、自分の作った my_dgemm と最適化 BLAS の dgemm_の速度を比べてみたが、その差は 1.5 倍ほど最適化 BLAS が速い、というものであった。 my_dgemm では一部にアセンブラを使っているが、これは自分の使った処理系がうまく命令をスケジューリングしていなかったからである。従って、プロセッサに特化されたような処理系であれば、アセンブラを使わずとも充分な最適化がなされることが予想される。つまり、良い処理系さえあれば、 C/C++言語だけで最適化 BLAS に更に近づくことが可能であろう。このような結果から、 ATLAS のように高級言語だけで最適化 BLAS に匹敵する BLAS を生成するという試みが成功していることも、納得できる。

LAPACK は最適化 BLAS と組み合わせることによって、非常に強力な数値計算ライブラリとなる。しかし、 FORTRAN のソースを f2c で変換したものをベースとしている CLAPACK のインターフェイスは、使いやすいものとは言えないであろう。まだ完全ではないようであるが、 LAPACK++ は C++で書かれており、クラスという形でインターフェイスを提供している。実際、 LAPACK 以外でも TNT(Template Numerical Toolkit) など、数値計算ライブラリにオブジェクト指向、或いはジェネリックプログラミングのインターフェイスを付ける流れはあるようである。

A 実行結果

A.1 dgemm
Testing double matrix multiplication.
alpha = 0.040329, beta = 0.0428009
N(a) : 1082x1025, N(b) : 1025x1071, c : 1082x1071
my_gemm : 6640msec

MKL GEMM : 4290msec
ratio : 1.54779
maximum error : 0
alpha = 0.231628, beta = 0.216263
N(a) : 1088x1054, T(b) : 1054x1034, c : 1088x1034
my_gemm : 6480msec
MKL GEMM : 5000msec
ratio : 1.296
maximum error : 0
alpha = 0.312836, beta = 0.138657
N(a) : 1059x1087, C(b) : 1087x1055, c : 1059x1055
my_gemm : 6760msec
MKL GEMM : 4830msec
ratio : 1.39959
maximum error : 0
alpha = 0.120361, beta = 0.191986
T(a) : 1016x1004, N(b) : 1004x1095, c : 1016x1095
my_gemm : 5880msec
MKL GEMM : 3950msec
ratio : 1.48861
maximum error : 0
alpha = 0.406311, beta = 0.189728
T(a) : 1009x1085, T(b) : 1085x1056, c : 1009x1056
my_gemm : 6370msec
MKL GEMM : 4510msec
ratio : 1.41242
maximum error : 0
alpha = 0.173157, beta = 0.0667572
T(a) : 1038x1046, C(b) : 1046x1036, c : 1038x1036
my_gemm : 6260msec
MKL GEMM : 4780msec
ratio : 1.30962
maximum error : 0
alpha = 0.113327, beta = 0.489944
C(a) : 1093x1078, N(b) : 1078x1009, c : 1093x1009
my_gemm : 6370msec
MKL GEMM : 4170msec
ratio : 1.52758
maximum error : 0
alpha = 0.436462, beta = 0.0402527
C(a) : 1011x1092, T(b) : 1092x1000, c : 1011x1000
my_gemm : 6100msec
MKL GEMM : 4340msec
ratio : 1.40553
maximum error : 0
alpha = 0.235825, beta = 0.0622101
C(a) : 1099x1035, C(b) : 1035x1066, c : 1099x1066
my_gemm : 6700msec
MKL GEMM : 4720msec
ratio : 1.41949

maximum error : 0

Testing float matrix multiplication.
alpha = 0.496094, beta = 0.761719
N(a) : 1096x1080, N(b) : 1080x1082, c : 1096x1082
my_gemm : 6980msec
MKL GEMM : 1970msec
ratio : 3.54315
maximum error : 6.10352e-005
alpha = 0.171875, beta = 0.917969
N(a) : 1017x1086, T(b) : 1086x1027, c : 1017x1027
my_gemm : 6150msec
MKL GEMM : 1870msec
ratio : 3.28877
maximum error : 1.90735e-005
alpha = 0, beta = 0.699219
N(a) : 1094x1022, C(b) : 1022x1031, c : 1094x1031
my_gemm : 0msec
MKL GEMM : 1920msec
ratio : 0
maximum error : 0
alpha = 0.902344, beta = 0.804688
T(a) : 1033x1008, N(b) : 1008x1040, c : 1033x1040
my_gemm : 5880msec
MKL GEMM : 1650msec
ratio : 3.56364
maximum error : 9.15527e-005
alpha = 0.167969, beta = 0.324219
T(a) : 1032x1082, T(b) : 1082x1004, c : 1032x1004
my_gemm : 6100msec
MKL GEMM : 1870msec
ratio : 3.26203
maximum error : 1.90735e-005
alpha = 0.175781, beta = 0.953125
T(a) : 1096x1063, C(b) : 1063x1043, c : 1096x1043
my_gemm : 6540msec
MKL GEMM : 2080msec
ratio : 3.14423
maximum error : 2.28882e-005
alpha = 0.46875, beta = 0.285156
C(a) : 1026x1075, N(b) : 1075x1083, c : 1026x1083
my_gemm : 6530msec
MKL GEMM : 1820msec
ratio : 3.58791
maximum error : 6.10352e-005
alpha = 0.496094, beta = 0.867188
C(a) : 1021x1036, T(b) : 1036x1028, c : 1021x1028
my_gemm : 5870msec
MKL GEMM : 1870msec

ratio : 3.13904
maximum error : 4.57764e-005
alpha = 0.746094, beta = 0.644531
C(a) : 1092x1032, C(b) : 1032x1023, c : 1092x1023
my_gemm : 6270msec
MKL GEMM : 1970msec
ratio : 3.18274
maximum error : 7.62939e-005

A.2 dgesv

A.2.1 標準 BLAS 使用

a : 189 x 189
b : 189 x 500
my_dgesv : 1100msec
dgesv_: 710msec
maximum error : 0
a : 500 x 500
b : 500 x 210
my_dgesv : 4890msec
dgesv_: 4890msec
maximum error : 0
a : 250 x 250
b : 250 x 250
my_dgesv : 930msec
dgesv_: 880msec
maximum error : 0
a : 20 x 20
b : 20 x 20
my_dgesv : 0msec
dgesv_: 0msec
maximum error : 0

A.2.2 最適化 BLAS 使用

a : 189 x 189
b : 189 x 500
my_dgesv : 280msec
dgesv_: 1040msec
maximum error : 1.63709e-011
a : 500 x 500
b : 500 x 210
my_dgesv : 610msec
dgesv_: 4610msec

maximum error : 1.64619e-010
a : 250 x 250
b : 250 x 250
my_dgesv : 110msec
dgesv_: 930msec
maximum error : 1.73372e-011
a : 20 x 20
b : 20 x 20
my_dgesv : 0msec
dgesv_: 0msec
maximum error : 1.13687e-013

B ソースコード

B.1 dgemm
template< class Scalar, class RowDim, class ColDim, class BlockDim > class block_t;
template< class Scalar, class BlockDim > class blockarray_t;
template< class Scalar, class Trans > class matrix_t;
///////////////////////////////////////////////////////////////////////////////
template< bool B > class Bool{};
typedef Bool< true >    True;
typedef Bool< false >   False;
template< size_t Dim >
class static_dim
{
public:
        typedef True    is_static;
        enum { dim = Dim };
        void    init( const static_dim< Dim > & ){}
};
class dynamic_dim
{
public:
        typedef False   is_static;
        const size_t    dim;
        dynamic_dim( void ) : dim( 0 ){}
        dynamic_dim( size_t d ) : dim( d ){}
        void    init( const dynamic_dim &d ){ const_cast< size_t & >( dim ) = d.dim; }
};
template< class Scalar, class RowDim, class ColDim, class BlockDim >
class block_t
{
public:
        void    init( RowDim r, ColDim c )
        {
                row.init( r );
                col.init( c );
                for ( size_t i = 0 ; i < row.dim ; ++i )
                        for ( size_t j = 0 ; j < col.dim ; ++j )
                                elements[i][j] = 0.0;

        }
        template< class Trans >
        void    init( RowDim ro, ColDim co, const matrix_t< Scalar, Trans > &m, size_t br1, size_t bc1 )
        {
                row.init( ro );
                col.init( co );
                for ( size_t i = 0, r = br1 * Dim * m.row_stride.dim ; i < row.dim ; ++i, r += m.row_stride.dim )
                        for ( size_t j = 0, c = bc1 * Dim * m.col_stride.dim ; j < col.dim ; ++j, c += m.col_stride.dim )
                                elements[i][j] = m.top[ r + c ];
        }
        enum { Dim = BlockDim::dim };
        RowDim  row;
        ColDim  col;
        Scalar  elements[Dim][Dim];
};
///////////////////////////////////////////////////////////////////////////////
class nontransposed
{
public:
        nontransposed( size_t ofs ) : row_stride( ofs ){}
        dynamic_dim     row_stride;
        static_dim< 1 > col_stride;
};
class transposed
{
public:
        transposed( size_t ofs ) : col_stride( ofs ){}
        static_dim< 1 > row_stride;
        dynamic_dim     col_stride;
};
template< class Scalar, class Trans >
class matrix_t : public Trans
{
public:
        Scalar  *const top;
        const size_t    row, col;
        matrix_t( Scalar *t, size_t r, size_t c, size_t ofs ) : Trans( ofs ), top( t ), row( r ), col( c ){}
        template< class RowDim, class ColDim, class BlockDim >
        void    assign( Scalar alpha, const block_t< Scalar, RowDim, ColDim, BlockDim > &blk, Scalar beta, size_t br1, size_t bc1 )
        {
                for ( size_t i = 0, r = br1 * blk.Dim * row_stride.dim ; i < blk.row.dim ; ++i, r += row_stride.dim )
                        for ( size_t j = 0, c = bc1 * blk.Dim * col_stride.dim ; j < blk.col.dim ; ++j, c += col_stride.dim )
                                top[ r + c ] = alpha * blk.elements[i][j] + beta * top[ r + c ];
        }
};
///////////////////////////////////////////////////////////////////////////////
#define LOOP2( x )      x x
#define LOOP4( x )      x x x x
#define LOOP5( x )      x x x x x
#define LOOP8( x )      LOOP2( LOOP4( x ) )
#define LOOP10( x )     LOOP2( LOOP5( x ) )
#define LOOP20( x )     LOOP2( LOOP10( x ) )
#define LOOP25( x )     LOOP5( LOOP5( x ) )

#define LOOP40( x )     LOOP4( LOOP10( x ) )
#define LOOP50( x )     LOOP5( LOOP10( x ) )
#define LOOP80( x )     LOOP8( LOOP10( x ) )
#define LOOP100( x )    LOOP10( LOOP10( x ) )
#ifdef _MSC_VER
void    __fastcall mul_block_asm( double *blk_c, const double *blk_a, const double *blk_b );
#endif
///////////////////////////////////////////////////////////////////////////////
template< class Scalar, class RowDimC, class ColDimC, class ColDimAB, class BlockDim >
static inline void      add_mul_sub( block_t< Scalar, RowDimC, ColDimC, BlockDim > &blk_c, const block_t< Scalar, RowDimC, ColDimAB, BlockDim > &blk_a, const block_t< Scalar, ColDimC, ColDimAB, BlockDim > &blk_b, False )
{
        for ( size_t i = 0 ; i < blk_a.row.dim ; ++i )
                for ( size_t j = 0 ; j < blk_b.row.dim ; ++j )
                        for ( size_t k = 0 ; k < blk_a.col.dim ; ++k )
                                blk_c.elements[i][j] += blk_a.elements[i][k] * blk_b.elements[j][k];
}
template< class Scalar, class RowDimC, class ColDimC, class ColDimAB, class BlockDim >
static inline void      add_mul_sub( block_t< Scalar, RowDimC, ColDimC, BlockDim > &blk_c, const block_t< Scalar, RowDimC, ColDimAB, BlockDim > &blk_a, const block_t< Scalar, ColDimC, ColDimAB, BlockDim > &blk_b, True )
{
        assert( ColDimAB::dim == 50 );
        for ( size_t i = 0 ; i < blk_a.row.dim ; ++i )
        {
                for ( size_t j = 0 ; j < blk_b.row.dim ; ++j )
                {
                        size_t  k = 0;
                        LOOP50
                        (
                                blk_c.elements[i][j] += blk_a.elements[i][k] * blk_b.elements[j][k];    ++k;
                        )
                }
        }
}
#ifdef _MSC_VER
typedef static_dim< 50 >        DefaultBlockDim;
template<>
static inline void      add_mul_sub< double, DefaultBlockDim, DefaultBlockDim, DefaultBlockDim, DefaultBlockDim >( block_t< double, DefaultBlockDim, DefaultBlockDim, DefaultBlockDim > &blk_c, const block_t< double, DefaultBlockDim, DefaultBlockDim, DefaultBlockDim > &blk_a, const block_t< double, DefaultBlockDim, DefaultBlockDim, DefaultBlockDim > &blk_b, True )
{
        mul_block_asm( &blk_c.elements[0][0], &blk_a.elements[0][0], &blk_b.elements[0][0] );
}
#endif
///////////////////////////////////////////////////////////////////////////////
template< class Scalar, class RowDimC, class ColDimC, class ColDimAB, class BlockDim >
static inline void      add_mul( block_t< Scalar, RowDimC, ColDimC, BlockDim > &blk_c, const block_t< Scalar, RowDimC, ColDimAB, BlockDim > &blk_a, const block_t< Scalar, ColDimC, ColDimAB, BlockDim > &blk_b )
{
        assert( blk_c.row.dim == blk_a.row.dim && blk_c.col.dim == blk_b.row.dim && blk_a.col.dim == blk_b.col.dim );
        add_mul_sub( blk_c, blk_a, blk_b, ColDimAB::is_static() );
}
///////////////////////////////////////////////////////////////////////////////
template< class Scalar, class RowDim, class ColDim, class BlockDim, class TransA, class TransB >
static inline void      my_gemm_block( block_t< Scalar, RowDim, ColDim, BlockDim > &blk_c, const matrix_t< Scalar, TransA > &a, const matrix_t< Scalar, TransB > &b, const size_t i, const size_t j )
{
        const size_t    bn_col = a.col / BlockDim::dim;
        const size_t    br_col = a.col % BlockDim::dim;

for ( size_t k = 0 ; k < bn_col ; ++k )
{
        static block_t< Scalar, RowDim, BlockDim, BlockDim >    blk_a;
        static block_t< Scalar, ColDim, BlockDim, BlockDim >    blk_b;
        blk_a.init( blk_c.row, BlockDim(), a, i, k );
        blk_b.init( blk_c.col, BlockDim(), b, j, k );
        add_mul( blk_c, blk_a, blk_b );
}
if ( br_col )
{
        static block_t< Scalar, RowDim, dynamic_dim, BlockDim > blk_a;
        static block_t< Scalar, ColDim, dynamic_dim, BlockDim > blk_b;
        blk_a.init( blk_c.row, br_col, a, i, k );
        blk_b.init( blk_c.col, br_col, b, j, k );
        add_mul( blk_c, blk_a, blk_b );
}
}
///////////////////////////////////////////////////////////////////////////////
template< class Scalar, class TransA, class TransB >
static void     my_gemm_sub( const Scalar alpha, const matrix_t< Scalar, TransA > &a, const matrix_t< Scalar, TransB > &b, const Scalar beta, matrix_t< Scalar, nontransposed > &c )
{
assert( c.get_row() == a.get_row() && c.get_col() == b.get_row() && a.get_col() == b.get_col() );
typedef static_dim< 50 >        BlockDim;
typedef block_t< Scalar, BlockDim, BlockDim, BlockDim >         center_block_t;
typedef block_t< Scalar, BlockDim, dynamic_dim, BlockDim >      right_block_t;
typedef block_t< Scalar, dynamic_dim, BlockDim, BlockDim >      bottom_block_t;
typedef block_t< Scalar, dynamic_dim, dynamic_dim, BlockDim >   rightbottom_block_t;
const size_t    bn_row_a = a.row / BlockDim::dim, bn_row_b = b.row / BlockDim::dim;
const size_t    br_row_a = a.row % BlockDim::dim, br_row_b = b.row % BlockDim::dim;
for ( size_t i = 0 ; i < bn_row_a ; ++i )
{
        for ( size_t j = 0 ; j < bn_row_b ; ++j )
        {
                static center_block_t   blk_c;
                blk_c.init( BlockDim(), BlockDim() );
                my_gemm_block( blk_c, a, b, i, j );
                c.assign( alpha, blk_c, beta, i, j );
        }
}
if ( br_row_b )
{
        const size_t    j = bn_row_b;
        for ( size_t i = 0 ; i < bn_row_a ; ++i )
        {
                static right_block_t    blk_c;
                blk_c.init( BlockDim(), br_row_b );
                my_gemm_block( blk_c, a, b, i, j );
                c.assign( alpha, blk_c, beta, i, j );
        }
}
if ( br_row_a )
{

const size_t    i = bn_row_a;
for ( size_t j = 0 ; j < bn_row_b ; ++j )
{
        static bottom_block_t   blk_c;
        blk_c.init( br_row_a, BlockDim() );
        my_gemm_block( blk_c, a, b, i, j );
        c.assign( alpha, blk_c, beta, i, j );
}
}
if ( br_row_a && br_row_b )
{
const size_t    i = bn_row_a;
const size_t    j = bn_row_b;
static rightbottom_block_t      blk_c;
blk_c.init( br_row_a, br_row_b );
my_gemm_block( blk_c, a, b, i, j );
c.assign( alpha, blk_c, beta, i, j );
}
}
template< class Scalar >
int     my_gemm( char *transa, char *transb, integer *m, integer *n, integer *k, Scalar *alpha, Scalar *a, integer *lda, Scalar *b, integer *ldb, Scalar *beta, Scalar *c, integer *ldc )
{
const char      None = 'N', Transpose = 'T', Conjugate = 'C', ta = toupper( *transa ), tb = toupper( *transb );
const bool      fTransA = ( ta != None ), fTransB = ( tb != None );
integer info = 0, nrowa, nrowb;
nrowa = ( fTransA ? *k : *m );
nrowb = ( fTransB ? *n : *k );
if ( ta != None && ta != Transpose && ta != Conjugate )
info = 1;
else if ( tb != None && tb != Transpose && tb != Conjugate )
info = 2;
else if ( *m < 0 )
info = 3;
else if ( *n < 0 )
info = 4;
else if ( *k < 0 )
info = 5;
else if ( *lda < max( 1, nrowa ) )
info = 8;
else if ( *ldb < max( 1, nrowb ) )
info = 10;
else if ( *ldc < max( 1, *m ) )
info = 13;
if ( info != 0 )
{
xerbla_( "MY_DGEMM ", &info );
return 0;
}
if ( *m == 0 || *n == 0 || ( ( *alpha == 0.0 || *k == 0 ) && *beta == 1.0 ) )
return 0;
if ( *alpha == 0.0 )
{
if ( *beta == 0.0 )
{
        for ( int j = 0 ; j < *n ; ++j )
                for ( int i = 0 ; i < *m ; ++i )
                        c[ j * *ldc + i ] = 0.0;
}

else
{
        for ( int j = 0 ; j < *n ; ++j )
                for ( int i = 0 ; i < *m ; ++i )
                        c[ j * *ldc + i ] *= *beta;
}
return 0;
}
if ( !fTransA )
{
if ( !fTransB )
        my_gemm_sub( *alpha, matrix_t< Scalar, nontransposed >( b, *n, *k, *ldb ), matrix_t< Scalar, transposed >( a, *m, *k, *lda ), *beta, matrix_t< Scalar, nontransposed >( c, *n, *m, *ldc ) );
else
        my_gemm_sub( *alpha, matrix_t< Scalar, transposed >( b, *n, *k, *ldb ), matrix_t< Scalar, transposed >( a, *m, *k, *lda ), *beta, matrix_t< Scalar, nontransposed >( c, *n, *m, *ldc ) );
}
else
{
if ( !fTransB )
        my_gemm_sub( *alpha, matrix_t< Scalar, nontransposed >( b, *n, *k, *ldb ), matrix_t< Scalar, nontransposed >( a, *m, *k, *lda ), *beta, matrix_t< Scalar, nontransposed >( c, *n, *m, *ldc ) );
else
        my_gemm_sub( *alpha, matrix_t< Scalar, transposed >( b, *n, *k, *ldb ), matrix_t< Scalar, nontransposed >( a, *m, *k, *lda ), *beta, matrix_t< Scalar, nontransposed >( c, *n, *m, *ldc ) );
}
return 0;
}
template int    my_gemm< real >( char *transa, char *transb, integer *m, integer *n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta, real *c, integer *ldc );
template int    my_gemm< doublereal >( char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc );

int     my_sgemm( char *transa, char *transb, integer *m, integer *n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta, real *c, integer *ldc )
{
return my_gemm( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc );
}
int     my_dgemm( char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc )
{
return my_gemm( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc );
}
///////////////////////////////////////////////////////////////////////////////
#ifdef _MSC_VER
__declspec( naked ) void        __fastcall mul_block_asm( double *blk_c, const double *blk_a, const double *blk_b )
{
// edx          : blk_a : call by register
// ebx, eax     : blk_b
// ecx          : blk_c : call by register
__asm
{
push    ebx
mov     ebx, DWORD PTR [esp+8]
push    esi
push    edi
add     edx, 16                                 ; 00000010H
mov     edi, 50                                 ; 00000032H
sub     ecx, 8
ROW:
mov     eax, ebx
mov     esi, 50                                 ; 00000032H

COL:
fld     QWORD PTR [eax+16+376]
fmul    QWORD PTR [edx+376]
fld     QWORD PTR [eax+16+368]
fmul    QWORD PTR [edx+368]
fld     QWORD PTR [eax+16+360]
fmul    QWORD PTR [edx+360]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+352]
fmul    QWORD PTR [edx+352]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+344]
fmul    QWORD PTR [edx+344]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+336]
fmul    QWORD PTR [edx+336]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+328]
fmul    QWORD PTR [edx+328]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+320]
fmul    QWORD PTR [edx+320]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+312]
fmul    QWORD PTR [edx+312]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+304]
fmul    QWORD PTR [edx+304]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+296]
fmul    QWORD PTR [edx+296]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+288]
fmul    QWORD PTR [edx+288]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+280]
fmul    QWORD PTR [edx+280]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+272]
fmul    QWORD PTR [edx+272]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+264]
fmul    QWORD PTR [edx+264]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+256]
fmul    QWORD PTR [edx+256]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+248]
fmul    QWORD PTR [edx+248]
faddp   ST(2), ST(0)

fld     QWORD PTR [eax+16+240]
fmul    QWORD PTR [edx+240]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+232]
fmul    QWORD PTR [edx+232]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+224]
fmul    QWORD PTR [edx+224]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+216]
fmul    QWORD PTR [edx+216]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+208]
fmul    QWORD PTR [edx+208]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+200]
fmul    QWORD PTR [edx+200]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+192]
fmul    QWORD PTR [edx+192]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+184]
fmul    QWORD PTR [edx+184]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+176]
fmul    QWORD PTR [edx+176]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16+168]
fmul    QWORD PTR [edx+168]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16+160]
fmul    QWORD PTR [edx+160]
add     eax, 400                                ; 00000190H
add     ecx, 8
dec     esi
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-248]
fmul    QWORD PTR [edx+152]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-256]
fmul    QWORD PTR [edx+144]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-264]
fmul    QWORD PTR [edx+136]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-272]
fmul    QWORD PTR [edx+128]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-280]

fmul    QWORD PTR [edx+120]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-288]
fmul    QWORD PTR [edx+112]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-296]
fmul    QWORD PTR [edx+104]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-304]
fmul    QWORD PTR [edx+96]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-312]
fmul    QWORD PTR [edx+88]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-320]
fmul    QWORD PTR [edx+80]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-328]
fmul    QWORD PTR [edx+72]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-336]
fmul    QWORD PTR [edx+64]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-344]
fmul    QWORD PTR [edx+56]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-352]
fmul    QWORD PTR [edx+48]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-360]
fmul    QWORD PTR [edx+40]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-368]
fmul    QWORD PTR [edx+32]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-376]
fmul    QWORD PTR [edx+24]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-416]
fmul    QWORD PTR [edx-16]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-384]
fmul    QWORD PTR [edx+16]
faddp   ST(2), ST(0)
fld     QWORD PTR [eax+16-408]
fmul    QWORD PTR [edx-8]
faddp   ST(1), ST(0)
fld     QWORD PTR [eax+16-392]
fmul    QWORD PTR [edx+8]
faddp   ST(2), ST(0)

                fld     QWORD PTR [edx]
                fmul    QWORD PTR [eax+16-400]
                faddp   ST(1), ST(0)
                fadd    QWORD PTR [ecx]
                faddp   ST(1), ST(0)
                fstp    QWORD PTR [ecx]
                jne     COL
                add     edx, 400                                ; 00000190H
                dec     edi
                jne     ROW
                pop     edi
                pop     esi
                pop     ebx
                ret     4
        }
}
#endif
///////////////////////////////////////////////////////////////////////////////
using std::cout;
using std::endl;
///////////////////////////////////////////////////////////////////////////////
void    gemm_( char *transa, char *transb, integer *m, integer *n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta, real *c, integer *ldc )
{
        SGEMM( transa, transb, ( int * )m, ( int * )n, ( int * )k, alpha, a, ( int * )lda, b, ( int * )ldb, beta, c, ( int * )ldc );
}
void    gemm_( char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc )
{
        DGEMM( transa, transb, ( int * )m, ( int * )n, ( int * )k, alpha, a, ( int * )lda, b, ( int * )ldb, beta, c, ( int * )ldc );
}
///////////////////////////////////////////////////////////////////////////////
template< class T > T   random( void ){ return random( ( T * )NULL ); }
double  random( double * )
{
        return ( ( rand() & 0xffff ) / ( double )0x10000 );
}
float   random( float * )
{
        return ( ( rand() & 0xff ) / ( float )0x100 );
}
///////////////////////////////////////////////////////////////////////////////
template< class Scalar >
Scalar  test( void )
{
        const size_t    Dim = 1100;
        const char      trans[] = "NTC";
        cout << "Testing " << typeid( Scalar ).name() << " matrix multiplication. " << endl << endl;
        for ( int ta = 0 ; trans[ta] ; ++ta )
        {
                for ( int tb = 0 ; trans[tb] ; ++tb )

{
static Scalar   a[ Dim * Dim ], b[ Dim * Dim ], c1[ Dim * Dim ], c2[ Dim * Dim ];
Scalar  alpha, beta;
integer m, n, k, lda, ldb, ldc;
char    transa, transb;
int     i;
transa = trans[ta];
transb = trans[tb];
alpha = random< Scalar >();
beta = random< Scalar >();
m = rand() % 100 + 1000;
n = rand() % 100 + 1000;
k = rand() % 100 + 1000;
lda = ldb = ldc = Dim;
cout    << "alpha = " << alpha << ", beta = " << beta << endl
        << transa << "(a) : " << m << "x" << k << ", "
        << transb << "(b) : " << k << "x" << n << ", "
        << "c : " << m << "x" << n << endl;

for ( i = 0 ; i < Dim * Dim ; ++i )
        a[i] = random< Scalar >();
for ( i = 0 ; i < Dim * Dim ; ++i )
        b[i] = random< Scalar >();
for ( i = 0 ; i < Dim * Dim ; ++i )
        c1[i] = c2[i] = random< Scalar >();

clock_t clk0, clk1, clkMy, clkMKL;
clk0 = clock();
my_gemm( &transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c1, &ldc );
clk1 = clock();
clkMy = clk1 - clk0;
cout << "my_gemm : " << ( double )clkMy / CLOCKS_PER_SEC * 1000.0 << "msec" << endl;
clk0 = clock();
gemm_( &transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c2, &ldc );
clk1 = clock();
clkMKL = clk1 - clk0;
cout << "MKL GEMM : " << ( double )clkMKL / CLOCKS_PER_SEC * 1000.0 << "msec" << endl;
cout << "ratio : " << ( double )clkMy / clkMKL << endl;

double  err = 0.0, tmp;
for ( i = 0 ; i < Dim * Dim ; ++i )
{
        tmp = fabs( c1[i] - c2[i] );
        err = ( tmp > err ? tmp : err );
}
cout << "maximum error : " << err << endl;
cout << endl;
}
}
cout << endl;
return Scalar();
}

int     main( void )
{
        srand( time( NULL ) );
        test< doublereal >();
        test< real >();
        return 0;
}
///////////////////////////////////////////////////////////////////////////////

B.2 dgesv
#define MKL
///////////////////////////////////////////////////////////////////////////////
#ifdef MKL
#include <mkl_blas.h>
#define idamax_( n, dx, incx )                                                  IDAMAX( ( int * )(n), dx, ( int * )(incx) )
#define dswap_( n, dx, incx, dy, incy )                                         DSWAP( ( int * )(n), dx, ( int * )(incx), dy, ( int * )(incy) )
#define dscal_( n, da, dx, incx )                                               DSCAL( ( int * )(n), da, dx, ( int * )(incx) )
#define dger_( m, n, alpha, x, incx, y, incy, a, lda )                          DGER( ( int * )(m), ( int * )(n), alpha, x, ( int * )(incx), y, ( int * )(incy), a, ( int * )(lda) )
#define dtrsm_( side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb )         DTRSM( side, uplo, transa, diag, ( int * )(m), ( int * )(n), alpha, a, ( int * )(lda), b, ( int * )(ldb) )
#define dgemm_( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )  DGEMM( transa, transb, ( int * )(m), ( int * )(n), ( int * )(k), alpha, a, ( int * )(lda), b, ( int * )(ldb), beta, c, ( int * )ldc )
#else
extern "C"
{
#include <fblaswr.h>
}
#endif
///////////////////////////////////////////////////////////////////////////////
#define A( row, col )   ( a[ (row) - 1 + ( (col) - 1 ) * *lda ] )
#define IPIV( x )       ( ipiv[ (x) - 1 ] )
///////////////////////////////////////////////////////////////////////////////
int     my_xerbla( char *srname, integer *info )
{
        printf( "** On entry to %6s, parameter number %2i had an illegal value\n", srname, *info );
        return 0;
}
///////////////////////////////////////////////////////////////////////////////
int     my_dlaswp( integer *n, doublereal *a, integer *lda, integer *k1, integer *k2, integer *ipiv, integer *incx )
{
        assert( *incx == 1 );
        for ( integer i = *k1 ; i <= *k2 ; ++i )
                if ( i != IPIV( i ) )
                        dswap_( n, &A( i, 1 ), lda, &A( IPIV( i ), 1 ), lda );
        return 0;
}
///////////////////////////////////////////////////////////////////////////////

int     my_dgetf2( integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info )
{
*info = 0;
if ( *m < 0 )
*info = -1;
else if ( *n < 0 )
*info = -2;
else if ( *lda < max( 1, *m ) )
*info = -4;
if ( *info != 0 )
{
integer tmp = -*info;
my_xerbla( "my_dgetf2", &tmp );
return 0;
}
if ( *m == 0 || *n == 0 )
return 0;
for ( integer j = 1 ; j <= min( *m, *n ) ; ++j )
{
integer tmp1, tmp2, tmp3;
integer jp;

// ピボットの選択

tmp1 = *m - j + 1;
tmp2 = 1;
jp = j - 1 + idamax_( &tmp1, &A( j, j ), &tmp2 );
IPIV( j ) = jp;
if ( A( jp, j ) != 0.0 )
{

// 第 j 行と第 jp 行の交換

        if ( jp != j )
                dswap_( n, &A( j, 1 ), lda, &A( jp, 1 ), lda );
        // (j+1,j) 成分から (*m,j) 成分を (j,j) 成分で割る
        if ( j < *m )
        {
                tmp1 = *m - j;
                doublereal      dtmp = 1.0 / A( j, j );
                tmp2 = 1;
                dscal_( &tmp1, &dtmp, &A( j + 1, j ), &tmp2 );
        }
}
else if ( *info == 0 )
{

// ピボット選択に失敗した場合

        *info = j;
}
if ( j < min( *m, *n ) )
{

// 列削除の分を引く

        tmp1 = *m - j;
        tmp2 = *n - j;
        doublereal      minus_one = -1.0;
        tmp3 = 1;
        dger_( &tmp1, &tmp2, &minus_one, &A( j + 1, j ), &tmp3, &A( j, j + 1 ), lda, &A( j + 1, j + 1 ), lda );
}
}
return 0;
}

///////////////////////////////////////////////////////////////////////////////
int     my_dgetrf( integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info )
{
*info = 0;
if ( *m < 0 )
*info = -1;
else if ( *n < 0 )
*info = -2;
else if ( *lda < max( 1, *m ) )
*info = -4;
if ( *info != 0 )
{
integer i = -*info;
my_xerbla( "my_dgetrf", &i );
return 0;
}
if ( *m == 0 || *n == 0 )
return 0;
const integer   NB = 50;
if ( NB <= 1 || NB >= min( *m, *n ) )
my_dgetf2( m, n, a, lda, ipiv, info );
else
{
for ( integer j = 1 ; j <= min( *m, *n ) ; j += NB )
{
integer iinfo, tmp1, tmp2, tmp3, tmp4;

// ブロックのサイズ

integer jb = min( min( *m, *n ) - j + 1, NB );
// (j,j)-(*m-1,j+jb-1) 行列を LU 分解
tmp1 = *m - j + 1;
my_dgetf2( &tmp1, &jb, &A( j, j ), lda, &IPIV( j ), &iinfo );
if ( *info == 0 && iinfo > 0 )
        *info = iinfo + j - 1;
for ( integer i = j ; i <= j + jb - 1 ; ++i )
        IPIV( i ) += j - 1;

// 1 列目から j-1 列目に、 j 行目から j+jb-1 行目の行交換を反映

tmp1 = j - 1;
tmp2 = j + jb - 1;
tmp3 = 1;
my_dlaswp( &tmp1, a, lda, &j, &tmp2, ipiv, &tmp3 );

// j+jb 列目以降が存在する場合

if ( j + jb <= *n )
{
        doublereal      one = 1.0;

// j+jb 列目から*n-1 列目に、 j 行目から j+jb-1 行目の行交換を反映

        tmp1 = *n - ( j + jb - 1 );
        tmp2 = j;
        tmp3 = j + jb - 1;
        tmp4 = 1;
        my_dlaswp( &tmp1, &A( 1, j + jb ), lda, &tmp2, &tmp3, ipiv, &tmp4 );
        // 上端部分 (j,j+jb)-(j+jb-1,*n-1) のブロックを計算
        tmp1 = *n - ( j + jb - 1 );
        dtrsm_( "Left", "Lower", "No transpose", "Unit", &jb, &tmp1, &one, &A( j, j ), lda, &A( j, j + jb ), lda );

// 右下部分が残っている場合

                        if ( j + jb <= *m )
                        {

// 列消去の分を引く

                                doublereal      minus_one = -1.0;
                                tmp1 = *m - ( j + jb - 1 );
                                tmp2 = *n - ( j + jb - 1 );
                                dgemm_( "No transpose", "No transpose", &tmp1, &tmp2, &jb, &minus_one, &A( j + jb, j ), lda, &A( j, j + jb ), lda, &one, &A( j + jb, j + jb ), lda );
                        }
                }
        }
}
return 0;
}
///////////////////////////////////////////////////////////////////////////////
int     my_dgetrs( char *trans, integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info )
{
assert( *trans == 'N' );
*info = 0;
if ( *n < 0 )
        *info = -2;
else if ( *nrhs < 0 )
        *info = -3;
else if ( *lda < max( 1, *n ) )
        *info = -5;
else if ( *ldb < max( 1, *n ) )
        *info = -8;
if ( *info != 0 )
{
        integer tmp = -*info;
        my_xerbla( "my_dgetrs", &tmp );
        return 0;
}
if ( *n == 0 || *nrhs == 0 )
        return 0;

integer tmp1, tmp2;
doublereal      one = 1.0;

// 右辺の行交換

tmp1 = 1;
tmp2 = 1;
my_dlaswp( nrhs, b, ldb, &tmp1, n, ipiv, &tmp2 );
// LY=B を解く
dtrsm_( "Left", "Lower", "No transpose", "Unit", n, nrhs, &one, a, lda, b, ldb );
// UX=Y を解く
dtrsm_( "Left", "Upper", "No transpose", "Non-unit", n, nrhs, &one, a, lda, b, ldb );
return 0;
}
///////////////////////////////////////////////////////////////////////////////
int     my_dgesv( integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info )
{
*info = 0;
if ( *n < 0 )

        *info = -1;
else if ( *nrhs < 0 )
        *info = -2;
else if ( *lda < max( 1, *n ) )
        *info = -4;
else if ( *ldb < max( 1, *n ) )
        *info = -7;
if ( *info != 0 )
{
        integer i = -*info;
        my_xerbla( "my_dgesv", &i );
        return 0;
}
my_dgetrf( n, n, a, lda, ipiv, info );
if ( *info == 0 )
        my_dgetrs( "No transpose", n, nrhs, a, lda, ipiv, b, ldb, info );
return 0;
}
///////////////////////////////////////////////////////////////////////////////
using std::cout;
using std::endl;
///////////////////////////////////////////////////////////////////////////////
void    test( integer n, integer nrhs )
{
const size_t    Dim = 500;
static doublereal       a1[ Dim * Dim ], b1[ Dim * Dim ], a2[ Dim * Dim ], b2[ Dim * Dim ];
static integer          ipiv[ Dim ];
integer         lda, ldb, info;
int             i;
clock_t         clk0, clk1;
lda = ldb = Dim;
for ( i = 0 ; i < Dim * Dim ; ++i )
{
        a1[i] = a2[i] = ( rand() % 0xffff ) / ( double )0x100;
        b1[i] = b2[i] = ( rand() % 0xffff ) / ( double )0x100;
}
cout    << "a : " << n << " x " << n << endl
        << "b : " << n << " x " << nrhs << endl << endl;
clk0 = clock();
my_dgesv( &n, &nrhs, a1, &lda, ipiv, b1, &ldb, &info );
clk1 = clock();
cout << "my_dgesv : " << ( clk1 - clk0 ) * 1000.0 / CLOCKS_PER_SEC << "msec" << endl;
clk0 = clock();
dgesv_( &n, &nrhs, a2, &lda, ipiv, b2, &ldb, &info );
clk1 = clock();
cout << "dgesv_: " << ( clk1 - clk0 ) * 1000.0 / CLOCKS_PER_SEC << "msec" << endl;
double  err = 0.0;
for ( i = 0 ; i < Dim * Dim ; ++i )
        err = max( err, max( fabs( a1[i] - a2[i] ), fabs( b1[i] - b2[i] ) ) );
cout << "maximum error : " << err << endl << endl;
}
///////////////////////////////////////////////////////////////////////////////

test( 189, 500 );
test( 500, 210 );
test( 250, 250 );
test( 20, 20 );
return 0;
}
///////////////////////////////////////////////////////////////////////////////