数値計算レポート

{G99P057-3 \atop 斎藤卓也}

出題日： 2001 年 6 月 6 日
提出期限： 2001 年 6 月 20 日
提出日： 2001 年 6 月 20 日

目次

1 問題 3 1.1 問題１ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 問題２ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 解答＆考察１ 3

3 実行結果１ 5

4 解答＆考察２ 6

5 実行結果２ 8

6 作成したプログラム１ 9 6.1 dgemm.c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 6.2 main.c dgemm.c チェック用 . . . . . . . . . . . . . . . . . . . . . . . . . . 18 6.3 MakeFile dgemm 用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

7 作成したプログラム 2 22 7.1 dgesv.c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 7.2 main.c dgesv.c チェック用 . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 7.3 Makefile dgesv 用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

1 問題

1.1 問題１

dgemm 関数を作れ（CLAPACK をインストールしながら仕様を見よ）。

1.2 問題２

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

dgesv (&n, &inc, A, &n piv, x, &n, &info)

を作れ。ただし、 BLAS を高速化したら高速になる形で実装せよ。つまり、標準化された BLAS の関数を使って dgesv (...) を作れ。

2 解答＆考察１

まず、私はインターネットから clapack.tgz をダウンロードし、 Vine Linux 2.1.5 上で LAPACK をインストールした。インストール自体は、 clapack.tgz を解凍したファイルで、 MakeFile を Linux 用と入れ替え、あとは Readme.Install の指示通りの作業をすることで完了した。

そして、 LAPACK の dgemm 関数の仕様を調べ、この関数と同じ動作をする dgemm 関数を作成した。本来なら、先日の行列高速化競争で作成したものを使用したかったのだが、この間作成したものは、 1000 × 1000 行列専用に徹底的に速くなるようにしたものだったので、どんなサイズの行列にもそのままでは対応できるようになっていない、また、 Vine Linux の gcc では SSE2 命令をアセンブラで使用した場合、そもそもコンパイルすら出来ない、今回のレポートでは特に高速化せよとの指示も無かったこと、そして、先週はレポートの数が大量にあり、今回のレポートに着手したのが遅く、４日程度しかなかったため、高速に演算できるよりも、確実に動くプログラムを作ろうと思い、そのようなプログラムを作成することにした。

作成した dgemm 関数の入力パラメータをチェックし、エラーを出力する部分のコードはそのまま使用した。これは、どちらにしても同じプログラムを組むことになるからである。

また、行列の演算を行う部分も、最初自分で組んでみようと思い、組み始めてから本物の dgemm 関数を見たところ、ほとんど同じようなプログラムになっていたので、改めて自分流に変更した。

といっても、やっていることは行列の掛け算なので、誰が組んでも似たようなプログラムになってしまうのだが、それは致し方ないところである。

また、一応高速化も図っている。具体的には、メモリの取り方を工夫し、前回と同じように、３２バイト領域にメモリ・アライメントを合わせて領域を確保している。

そして、行列の演算を少しでも高速化するために、掛け算の部分では、 A × B の B の行列を転地行列化している。

実際に速度はチェックしていないが、これで本来ある標準 BLAS のものよりは高速に演算できるであろう。しかし、もっとも高速化に寄与するブロック化や SSE2 は使用していないため、大した速度ではないと思われる。

また、以下に実行結果を示すが、ここでは

A = ( \begin{matrix}1 & 2 \\ 3 & 4 \end{matrix} )

B = ( \begin{matrix}5 & 6 \\ 7 & 8 \end{matrix} )

とし、

AB = ( \begin{matrix}19 & 22 \\ 43 & 50 \end{matrix} )

tAB = ( \begin{matrix}26 & 30 \\ 38 & 44 \end{matrix} )

A^{t}B = ( \begin{matrix}17 & 23 \\ 39 & 53 \end{matrix} )

tA^{t}B = ( \begin{matrix}23 & 31 \\ 34 & 46 \end{matrix} )

の計算結果を求めてみている。それぞれ、きちんと正しい結果が得られていることが確認できる。

3 実行結果１

A=
1.000000 2.000000
3.000000 4.000000
B=
5.000000 6.000000
7.000000 8.000000
A x B =
19.000000 22.000000
43.000000 50.000000
A' x B =
26.000000 30.000000
38.000000 44.000000
A x B' =
17.000000 23.000000
39.000000 53.000000
A' x B' =
23.000000 31.000000
34.000000 46.000000

4 解答＆考察２

次に、 dgesv 関数だが、こちらは標準 BLAS 関数を使って、 dgesv 関数を作れ、とのことだったので、まずはいろいろと LU 分解法で連立方程式を解くときに使用するのに良い標準 BLAS はどのようなものがあるのかを調べることからはじめた。

いろいろと調べた結果、 dtrsm という標準 BLAS 関数が丁度使用できると思い、この標準 BLAS 関数を使用した dgesv 関数を作成した。

この dtrsm という関数は、

op( A )*X = alpha*B,   or   X*op( A ) = alpha*B

という形の連立方程式を解いてくれる（ただし、 A は LU 分解されたものが入っている必要がある）。

もっと厳密に言うと、この A という行列は、 LU 分解で得られた、 L 行列と U 行列を１つにした行列で（そうでなくてもよい）、 A が L の成分を含んでいると考えるか（L の対角線上は全て１だと仮定することもできる）、あるいは、 U の成分を含んでいると考えるかを指定することにより、前進消去、後退代入の作業をどちらでもできるようになっているという、非常にすぐれものなのである。

元々、前進消去、後退代入の作業自体は大した作業ではないので、自分でプログラムを組んでもすぐに出来てしまうのだが、それよりもここでは BLAS 関数を使用することにより、 BLAS を高速なものに変えてやれば、その処理も高速化されるというところにポイントがある。

dgesv 関数の方も、 dgemm 関数と同様に、入力パラメータのエラーチェックに関しては、そのまま前の関数のものを使用することにした。これはあえて自分で作る理由もないからである。

基本的には、 LU 分解を行う部分を自分でプログラムしているのだが、前回の課題同様、特に難しいプログラムではない。

ただ、今回のプログラムでは、 dgesv 関数では、引数の行列 A が、実行が終わった後には LU 分解した

A = L 行列（対角線の１と右上部分は除く）＋ U 行列の形（左下の部分は除く）

の形で返さないといけないので、それに沿う形でプログラムを変更した。

また、このようにしておくと、 dtrsm 関数で LU 分解を使用して連立方程式を解かせるときにも、 uplo と diag パラメータを変更し、２度呼んでやるだけで前進消去、後退代入を行えるので、何かと都合が良い。

その他のポイントとしては、 Ax = b という方程式を解く場合、 dgesv 関数では、この b が複数列の行列の場合も対応できるようになっているので、それにも対応できるように変更してある。

この辺はただ、 b を１列から複数列対応にするだけなので、別に難しいところはなにもない。しかも、 dtrsm 関数も、そのまま対応しているので、非常に簡単にプログラムを組むことが出来た。

以上がプログラムの概要である。

次に、以下に実行結果を示してあるが、これは、

( \begin{matrix}2 & 6 & 4 \\ 3 & 7 & 2 \\ 1 & 2 & 3 \end{matrix} ) ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ) = ( \begin{matrix}8 \\ 6 \\ 7 \end{matrix} )

という連立方程式を解かせた結果である。実行結果の A=PLU=で示されているのが、 dgesv 関数を実行後に A 行列に代入された、 LU 分解後の値である。左下部分には L 行列（対角線上の１と右上除く）と U 行列（左下部分除く）を加えて一つにしたものである。

また、 IPIV はピボット選択で、どのように変更されたかが入っている。この例では、上から２行と１行が入れかえられており、３行はそのままになっているということを示している。

そして、この連立方程式を解くと、解は

( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ) = ( \begin{matrix}3 \\ −1 \\ 2 \end{matrix} )  (1)

であり、 dgesv 関数で解かせた結果は、もちろんこれと一致している。

5 実行結果２

A=
2.000000 6.000000 4.000000
3.000000 7.000000 2.000000
1.000000 2.000000 3.000000
b=
8.000000
6.000000
7.000000
A=PLU=
3.000000 7.000000 2.000000
0.666667 1.333333 2.666667
0.333333 -0.250000 3.000000
IPIV=
2
1
3
x=
3.000000
-1.000000
2.000000

6 作成したプログラム１

6.1 dgemm.c

/*-----------------------------------------------------------
dgemm_
2001-6-17
Programmed by Takuya Saito
-------------------------------------------------------------
Purpose
=======
DGEMM  performs one of the matrix-matrix operations
   C := alpha*op( A )*op( B ) + beta*C,
where  op( X ) is one of
   op( X ) = X   or   op( X ) = X',
alpha and beta are scalars, and A, B and C are matrices, with op( A )
an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
Parameters
==========
TRANSA - CHARACTER*1.
         On entry, TRANSA specifies the form of op( A ) to be used in
         the matrix multiplication as follows:
            TRANSA = 'N' or 'n',  op( A ) = A.
            TRANSA = 'T' or 't',  op( A ) = A'.
            TRANSA = 'C' or 'c',  op( A ) = A'.

Unchanged on exit.
TRANSB - CHARACTER*1.
On entry, TRANSB specifies the form of op( B ) to be used in
the matrix multiplication as follows:
   TRANSB = 'N' or 'n',  op( B ) = B.
   TRANSB = 'T' or 't',  op( B ) = B'.
   TRANSB = 'C' or 'c',  op( B ) = B'.
Unchanged on exit.
M      - INTEGER.
On entry,  M  specifies  the number  of rows  of the  matrix
op( A )  and of the  matrix  C.  M  must  be at least  zero.
Unchanged on exit.
N      - INTEGER.
On entry,  N  specifies the number  of columns of the matrix
op( B ) and the number of columns of the matrix C. N must be
at least zero.
Unchanged on exit.
K      - INTEGER.
On entry,  K  specifies  the number of columns of the matrix
op( A ) and the number of rows of the matrix op( B ). K must
be at least  zero.
Unchanged on exit.
ALPHA  - DOUBLE PRECISION.
On entry, ALPHA specifies the scalar alpha.

Unchanged on exit.
A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
part of the array  A  must contain the matrix  A,  otherwise
the leading  k by m  part of the array  A  must contain  the
matrix A.
Unchanged on exit.
LDA    - INTEGER.
On entry, LDA specifies the first dimension of A as declared
in the calling (sub) program. When  TRANSA = 'N' or 'n' then
LDA must be at least  max( 1, m ), otherwise  LDA must be at
least  max( 1, k ).
Unchanged on exit.
B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
part of the array  B  must contain the matrix  B,  otherwise
the leading  n by k  part of the array  B  must contain  the
matrix B.
Unchanged on exit.
LDB    - INTEGER.
On entry, LDB specifies the first dimension of B as declared
in the calling (sub) program. When  TRANSB = 'N' or 'n' then

LDB must be at least  max( 1, k ), otherwise  LDB must be at
least  max( 1, n ).
Unchanged on exit.
BETA   - DOUBLE PRECISION.
On entry,  BETA  specifies the scalar  beta.  When  BETA  is
supplied as zero then C need not be set on input.
Unchanged on exit.
C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
Before entry, the leading  m by n  part of the array  C must
contain the matrix  C,  except when  beta  is zero, in which
case C need not be set on entry.
On exit, the array  C  is overwritten by the  m by n  matrix
( alpha*op( A )*op( B ) + beta*C ).
LDC    - INTEGER.
On entry, LDC specifies the first dimension of C as declared
in  the  calling  (sub)  program.   LDC  must  be  at  least
max( 1, m ).
Unchanged on exit.
Level 3 Blas routine.
-- Written on 8-February-1989.
Jack Dongarra, Argonne National Laboratory.
Iain Duff, AERE Harwell.
Jeremy Du Croz, Numerical Algorithms Group Ltd.
Sven Hammarling, Numerical Algorithms Group Ltd.
Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not

   transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
   and  columns of  A  and the  number of  rows  of  B  respectively.
*/
#include "f2c.h"
#define A(i, j) a[(i) + (j)*( *lda)]
#define B(i, j) b[(i) + (j)*( *ldb)]
#define C(i, j) c[(i) + (j)*( *ldc)]
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)
{
extern logical lsame_(char*, char*);
extern int xerbla_(char*, integer*);
static logical nota, notb;
static integer info;
static integer ncola;
static integer nrowa, nrowb;
double **a2, **b2, **c2;
double *a2r, *b2r, *c2r;
double temp;
int gyouA, retuA, gyouB, retuB, i, j, w;
/*

エラー処理

*/
nota = lsame_(transa, "N");
notb = lsame_(transb, "N");
if (nota) {
    nrowa = *m;
    ncola = *k;

} else {
    nrowa = *k;
    ncola = *m;
}
if (notb) {
    nrowb = *k;
} else {
    nrowb = *n;
}
info = 0;
if (!nota && !lsame_(transa, "C") && !lsame_(transa, "T")) {
    info=1;
} else if (!notb && !lsame_(transb, "C") && !lsame_(transb, "T")) {
    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_("DGEMM ", &info);
    return 0;
}
/*

A と B の行数・列数をセット

*/
if (nota) {

    gyouA = *m;
    retuA = *k;
} else {
    gyouA = *k;
    retuA = *m;
}
if (notb) {
    gyouB = *k;
    retuB = *n;
} else {
    gyouB = *n;
    retuB = *k;
}
/*
  行列用メモリ領域確保 (a2, b2, c2)

３２バイト境界にアライメントを合わせてメモリを確保する。

*/
a2 = (double**) malloc ( sizeof (double*) * gyouA );
b2 = (double**) malloc ( sizeof (double*) * gyouB );
c2 = (double**) malloc ( sizeof (double*) * gyouA );
a2r = (double*) malloc (sizeof (double) * gyouA * retuA + 31);
b2r = (double*) malloc (sizeof (double) * gyouA * retuA + 31);
c2r = (double*) malloc (sizeof (double) * gyouA * retuA + 31);
a2[0] = (double*) (((unsigned)a2r + 31) & (-32));
b2[0] = (double*) (((unsigned)b2r + 31) & (-32));
c2[0] = (double*) (((unsigned)c2r + 31) & (-32));
for (i=1; i < gyouA; i++) {
    a2[i] = a2[i-1] + retuA;
    c2[i] = c2[i-1] + retuB;
}
for (i=1; i < gyouB; i++) {
    b2[i] = b2[i-1] + retuB;
}
/*

行列の要素をコピーする

*/
/* A */
if (nota) {
    for (i=0; i < gyouA; i++) {
        for (j=0; j < retuA; j++) {
            a2[i][j] = A(i, j);
        }
    }
} else {
    for (i=0; i < gyouA; i++) {
        for (j=0; j < retuA; j++) {
            a2[i][j] = A(j, i);
        }
    }
}
/* B */
if (notb) {
    for (i=0; i < gyouB; i++) {
        for (j=0; j < retuB; j++) {
            b2[j][i] = B(i, j);
        }
    }
} else {
    for (i=0; i < gyouB; i++) {
        for (j=0; j < retuB; j++) {
            b2[j][i] = B(j, i);
        }
    }
}
/*
  C := alpha * op(A)*op(B) + beta * C
*/
for (i=0; i < gyouA; i++) {
    for (j=0; j < retuB; j++) {
        temp = 0.0;
        for (w=0; w < retuA; w++) {

            temp += a2[i][w] * b2[j][w];
        }
        c2[i][j] = *alpha * temp + *beta * C(i, j);
    }
}
for (i=0; i < gyouA; i++) {
    for (j=0; j < retuB; j++) {
        C(i, j) = c2[i][j];
    }
}
/*

行列用メモリ領域解放

*/
free(a2r);
free(b2r);
free(c2r);
free(a2);
free(b2);
free(c2);
return 0;
}

6.2 main.c dgemm.c チェック用

/*
dgemm_test program
2001-6-20
Programmed by Takuya Saito
*/
#include <stdio.h>
#include <stdlib.h>
#include "f2c.h"
void disp_matrix(doublereal* a);
int main ()
{
double a[2*2];
double b[2*2];
double c[2*2];
int i, j;
char transa, transb;
integer m, n, k, lda, ldb, ldc;
doublereal alpha, beta;
a[0] = 1.; a[2] = 2.;
a[1] = 3.; a[3] = 4.;
b[0] = 5.; b[2] = 6.;
b[1] = 7.; b[3] = 8.;
for (i=0; i < 4; i++)
    c[i] = 1.0;
m = 2;
n = 2;
k = 2;
lda = 2;
ldb = 2;
ldc = 2;
alpha = 1.0;

beta = 0.0;
transa = 'n';
transb = 'n';
printf("A=\n");
disp_matrix(a);
printf("\nB=\n");
disp_matrix(b);
/* A * B */
dgemm_(&transa, &transb, &m, &n, &k, &alpha, a, &lda,
       b, &ldb, &beta, c, &ldc);
printf("\nA x B =\n");
disp_matrix(c);
/* A' * B */
transa = 't';
dgemm_(&transa, &transb, &m, &n, &k, &alpha, a, &lda,
       b, &ldb, &beta, c, &ldc);
printf("\nA' x B =\n");
disp_matrix(c);
/* A * B' */
transa = 'n';
transb = 't';
dgemm_(&transa, &transb, &m, &n, &k, &alpha, a, &lda,
       b, &ldb, &beta, c, &ldc);
printf("\nA x B' =\n");
disp_matrix(c);
/* A' * B' */
transa = 't';
transb = 't';
dgemm_(&transa, &transb, &m, &n, &k, &alpha, a, &lda,

       b, &ldb, &beta, c, &ldc);
printf("\nA' x B' =\n");
disp_matrix(c);
return 0;
}
/*

行列表示関数

*/
void disp_matrix(doublereal* a)
{
int i, j;
for (i=0; i < 2; i++) {
    for (j=0; j < 2; j++) {
        printf("%f ",a[i+2*j]);
    }
    printf("\n");
}
}

6.3 MakeFile dgemm 用

LIB=/usr/local/lib/lapack_LINUX.a /usr/local/lib/blas_LINUX.a /usr/local/lib/libF77.a
main: main.o dgemm.o $(LIB)
    cc -o main main.o dgemm.o $(LIB)
main.o:main.c
    cc -c main.c
dgemm.o: dgemm.c
    cc -c dgemm.c

7 作成したプログラム 2

7.1 dgesv.c

/*---------------------------------------------------
dgesv
2001-6-20
Programmed by Takuya Saito
-----------------------------------------------------

Purpose
=======
DGESV computes the solution to a real system of linear equations
   A * X = B,
where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
The LU decomposition with partial pivoting and row interchanges is
used to factor A as
   A = P * L * U,
where P is a permutation matrix, L is unit lower triangular, and U is
upper triangular.  The factored form of A is then used to solve the
system of equations A * X = B.
Arguments
=========
N       (input) INTEGER
        The number of linear equations, i.e., the order of the
        matrix A.  N >= 0.
NRHS    (input) INTEGER
        The number of right hand sides, i.e., the number of columns
        of the matrix B.  NRHS >= 0.
A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
        On entry, the N-by-N coefficient matrix A.

On exit, the factors L and U from the factorization
A = P*L*U; the unit diagonal elements of L are not stored.
LDA     (input) INTEGER
The leading dimension of the array A.  LDA >= max(1,N).
IPIV    (output) INTEGER array, dimension (N)
The pivot indices that define the permutation matrix P;
row i of the matrix was interchanged with row IPIV(i).
B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
On entry, the N-by-NRHS matrix of right hand side matrix B.
On exit, if INFO = 0, the N-by-NRHS solution matrix X.
LDB     (input) INTEGER
The leading dimension of the array B.  LDB >= max(1,N).
INFO    (output) INTEGER
= 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
      has been completed, but the factor U is exactly
      singular, so the solution could not be computed.
=====================================================================
*/
#include <stdio.h>
#include <math.h>
#include "f2c.h"
#define A(i, j) a[(i)-1 + ((j)-1) * (*lda)]
#define B(i, j) b[(i)-1 + ((j)-1) * (*lda)]
#define IPIV(I) ipiv[(I)-1]
int dgesv_(integer* n, integer* nrhs, doublereal* a, integer* lda,
integer* ipiv, doublereal* b, integer* ldb, integer* info)
{
extern int xerbla_(char*, integer*);
extern int dtrsm_(char*, char*, char*, char*, integer*, integer*,

                  doublereal*, doublereal*, integer*, doublereal*, integer*);
int i, j, k, w, q, q2, tempi;
double temp, alpha;
char side, uplo, transa, diag;
double A2[9];
*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) {
    i = -(*info);
    xerbla_("DGESV ", &i);
    return 0;
}
for (i=1; i <= *n; i++) {
    IPIV(i) = i;
}
/*
  LU 分解
*/
for (i=1; i < *n; i++) {
    /*

部分ピボット選択

    */
    q2 = i;
    for (q=i+1; q <= *n; q++) {
        if (fabs((double)A(q2, i)) < fabs((double)A(q, i))) {

    q2 = q;
}
}
if (A(q2, i) == 0.0) {
printf("Can not solve!\n");
return 0;
}
if (q2 != i) {
for (w=i; w <=*n ; w++) {
    temp = A(i, w);
    A(i, w) = A(q2, w);
    A(q2, w) = temp;
}
for (w=1; w <= *nrhs; w++) {
    temp = B(i, w);
    B(i, w) = B(q2, w);
    B(q2, w) = temp;
}
tempi = IPIV(i);
IPIV(i) = q2;
IPIV(q2) = tempi;
}

/*

LU 分解をする

*/
for (j=i+1; j <= *n; j++) {
/* j 行目から、 i 行目の a[j][i] 倍を引く。 */
temp = A(j, i) / A(i, i);
/*L[j-1][i-1] = temp;*/
for (k=i; k <= *n; k++) {
    A(j, k) -= temp * A(i, k);
}
A(j, i) = temp;
}
}

/*

標準ＢＬＡＳを使って連立方程式を解く

*/
side = 'L';
uplo = 'L';
transa = 'N';
diag = 'u';
alpha = 1.0;
dtrsm_(&side, &uplo, &transa, &diag, n, nrhs, &alpha, a, n, b, n);
uplo = 'U';
diag = 'n';
dtrsm_(&side, &uplo, &transa, &diag, n, nrhs, &alpha, a, n, b, n);
return 0;
}

7.2 main.c dgesv.c チェック用

/*

dgesv_チェック用プログラム

2001-6-20
Programmed by Takuya Saito
*/
#include <stdio.h>
#include <math.h>
#include "f2c.h"
void disp_matrix(doublereal* a, int n);
int main()
{
int i, j;
integer n, nrhs, lda, ldb, info;
integer ipiv[3];
doublereal a[3*3];
doublereal b[3];
a[0] = 2.0; a[3] = 6.0; a[6] = 4.0;
a[1] = 3.0; a[4] = 7.0; a[7] = 2.0;
a[2] = 1.0; a[5] = 2.0; a[8] = 3.0;
b[0] = 8.0; b[1] = 6.0; b[2] = 7.0;
n=3;
nrhs = 1;
lda = 3;
ldb = 3;
/* 行列 A 表示 */
printf("A=\n");
disp_matrix(a, 3);
/* b 表示 */
printf("\nb=\n");
for (i=0; i < 3; i++)

    printf("%f\n", b[i]);
dgesv_(&n,&nrhs, a, &lda, ipiv, b, &ldb, &info);

/* LU 分解後の行列 A 表示 */

printf("\nA=PLU=\n");
disp_matrix(a, 3);
/* IPIV 表示 */
printf("\nIPIV=\n");
for (i=0; i < 3; i++)
    printf("%d\n", ipiv[i]);
/* 答え表示 */
printf("\nx=\n");
for (i=0; i < 3; i++)
    printf("%f\n", b[i]);
return 0;
}
/*

行列表示関数

*/
void disp_matrix(doublereal* a, int n)
{
int i, j;
for (i=0; i < n; i++) {
    for (j=0; j < n; j++) {
        printf("%f ",a[i+n*j]);
    }
    printf("\n");
}
}

7.3 Makefile dgesv 用

LIB=/usr/local/lib/lapack_LINUX.a /usr/local/lib/blas_LINUX.a /usr/local/lib/libF77.a
main: main.o dgesv.o dtrsm.o $(LIB)
    cc -o main main.o dgesv.o dtrsm.o $(LIB)
main.o:main.c
    cc -c main.c
dgesv.o: dgesv.c
    cc -c dgesv.c
dtrsm.o: dtrsm.c
    cc -c dtrsm.c