C言語からLAPACKを呼んで逆行列を求める(LU分解する)
数値解析の授業で逆行列は求めちゃダメって言われたけど、気にしない気にしない。
実行列の逆行列を求める
Fortranの DGETRF でLU分解した後、 DGETRI で逆行列を求める。
/* * hoge.c * SIZE*SIZE型の実行列の逆行列を計算 * 元の行列は * (2 3 ) * (1 0.5) */ #define SIZE 2 #include <stdio.h> int main(void){ int m = SIZE ; // 行のサイズ int n = SIZE ; // 列のサイズ int lda = SIZE ; // mと同じ値 double A[SIZE*SIZE] ; // m x n の行列成分。この行列の逆行列を求める。 A[0] = 2.0 ;A[2] = 3.0; A[1] = 1.0 ;A[3] = 0.5; int info ; // 計算が成功すれば0を返す int ipiv[SIZE] ; // 要素数はm,nのうち小さい方とする int lwork = SIZE ; // nと同じ値 double work[SIZE] ; // 要素数はlworkと同じ値 // LAPACKのdgetrfサブルーチンを呼んで、行列AをLU分解 // 引数は全て参照渡し dgetrf_( &m, &n, A, &lda, ipiv, &info); // LU分解後の行列から逆行列を求める // 逆行列は元の配列Aに入る dgetri_( &n, A, &lda, ipiv, work, &lwork, &info); printf("%+8.5lf %+8.5lf\n", A[0], A[2]); printf("%+8.5lf %+8.5lf\n", A[1], A[3]); }
コンパイルと実行結果
$ gcc hoge.c -llapack -lblas -lm $ a.out -0.25000 +1.50000 +0.50000 -1.00000
引数の詳細については、下手な解説を見るより dgetrf と dgetri を読むのが一番わかりやすい。
CとFortranの配列形式の違いについては、 CとFortranで行列の添字が異なる点への注意喚起 に書いた。
複素行列の逆行列を求める
Fortarnの
ZGETRF
でLU分解した後、
ZGETRI
で逆行列を求める。
実行列の場合との違いは、行列Aと配列workの型がdouble
からdouble _Complex
に変わったところ。
/* * hoge.c * SIZE*SIZE型の複素行列の逆行列を計算 * 元の行列は * (0.5+i 1.0+0.5i) * (2.0 1.0 ) */ #define SIZE 2 // 2*2型の行列 #include <stdio.h> #include <complex.h> int main(void){ int m = SIZE ; // 行のサイズ int n = SIZE ; // 列のサイズ int lda = SIZE ; // mと同じ値 double _Complex A[SIZE*SIZE] ; // m x n の行列成分。この行列の逆行列を求める。 A[0]= 0.5+I ; A[2]= 1.0+0.5*I; A[1]= 2.0 ; A[3]= 1.0; int info ; // 計算が成功すれば0を返す int ipiv[SIZE] ; // 要素数はm,nのうち小さい方とする int lwork = SIZE ; // nと同じ値 double _Complex work[SIZE] ; // 要素数はlworkと同じ値 // LAPACKのzgetrfサブルーチンを呼んで、行列AをLU分解 // 引数は全て参照渡し zgetrf_( &m, &n, A, &lda, ipiv, &info); // LU分解後の行列から逆行列を求める // 逆行列は元の配列Aに入る zgetri_( &n, A, &lda, ipiv, work, &lwork, &info); printf("%+8.5lf%+8.5lfI %+8.5lf%+8.5lfI\n", creal(A[0]), cimag(A[0]), creal(A[2]), cimag(A[2])); printf("%+8.5lf%+8.5lfI %+8.5lf%+8.5lfI\n", creal(A[1]), cimag(A[1]), creal(A[3]), cimag(A[3])); }
コンパイルと実行結果
$ gcc hoge.c -llapack -lblas -lm $ a.out -0.66667-0.00000I +0.66667+0.33333I +1.33333+0.00000I -0.33333-0.66667I
引数の詳細は zgetrf と zgetri を読むべし。 またCとFortranの配列形式の違いについては、 CとFortranで行列の添字が異なる点への注意喚起 を読むべし。