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で行列の添字が異なる点への注意喚起 を読むべし。