Windows での LAPACK の使い方（Windows 上の Cygwin を使用）

LAPACK とは，行列に関する種々の問題（連立１次方程式，固有値問題，などなど多数）． FORTRAN で書かれています．

このページでは，LAPACK の下記の機能を C プログラムから呼び出す方法を，プログラム例を使って説明する．

dgeev()関数（固有値と固有ベクトル）

BLAS とは？

BLAS(Basic Linear Algebra Subprograms)とは，行列演算，ベクトル演算の機能をもったプログラム群です． GotoBLAS や ATLAS は，無料で使えるソフトウェアの中では，BLAS 実装の決定版といってよいでしょう．※ 他にも，優れた実装が複数あります．

＜BLAS の主な機能（ごく一部を紹介）＞

Level 1 ベクトルとベクトルの演算
- DOT : 内積
- AXPY : AXPY 演算 ( y <- ax + y の形など)
- NORM : ノルム　など
Level 2 行列とベクトルと計算
- 行列とベクトルの積 ( y <- Ax )
- 行列の rank-1 更新 ( A <- A + xy' )
Level 3 行列同士の演算
- 行列と行列の積 ( Z <- XY )

* GotoBLAS はプログラムで，しかも，無料で使える．但し，GotoBLAS はフリーソフトウェアではない．「再配布や，商用・業務利用は原則不可」と明記されている．再配布不可ですから，ソースコードやビルドしたバイナリを人にあげることはできません（パソコンを借りてインストールしてあげる，ということも難しい，という意味です）．

* GotoBLAS のライセンス条項に「不安」がある場合は，ATLAS を使うことも検討してください．（決して違法なことはしないように）．

LAPACK とは？

LAPACK とは，行列に関する種々の計算（連立１次方程式，固有値問題，などなど多数）のソフトウェア． LAPACK は，BLAS の機能を利用して，種々の機能を実現します．（BLAS を速いものに入れ替えれば，LAPACK も速くなるということです）．

前準備

前もってインストールしておくべきソフトウェア

BLAS 及び LAPACK のビルドとインストール
- GotoBLAS を使う場合
  「Windows で GotoBLAS と CBLAS をビルドとインストール」の Web ページに従って， GotoBLAS のインストールが終わっていること．さらに，
  「Windows で LAPACK バージョン 3.2.1 をビルドとインストール」の Web ページに従って， LAPACK のインストールが終わっていること．ライセンス条項を確認しておくこと．
- ATLAS を使う場合
  「Windows で ATLAS と LAPACK をビルドとインストール」の Web ページに従って， ATLAS と LAPACK のインストールが終わっていること．ライセンス条項を確認しておくこと．
GSL のビルドとインストール
「Windows に GSL (GNU Scientific Library) をインストール」の Web ページに従って， GSL (GNU Scientific Library) のインストールが終わっていること．

プログラムの見本

【要点】

C のプログラム
LAPACK 及び BLAS のオブジェクトファイルとリンクする
CLAPACK は使わない．
使用する LAPACK の機能
- dgeev()関数（固有値と固有ベクトル）

#include >stdio.h<

// CLAPACK ではなく LAPACK を直接呼び出す．

// CLAPACK については，次の URL を見よ．
//   see http://www.netlib.org/clapack/clapack.h

// ?geev : simple driver for eigenvalues/vectors

// 決まり文句 (CLAPACK との互換用）

#define doublereal double
#define integer long int

extern
  integer dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *a,
         integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
         integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
         integer *lwork, integer *info);

integer eigenvalues( integer n, doublereal *a, doublereal *wr, doublereal *wi ) {
  /* LAPACK の _dgeev() を使って固有値（だけ）を求める */

  integer n3 = n * n * n;
  integer info;

  doublereal *vl = (doublereal *)calloc(sizeof(doublereal), n * n);
  doublereal *vr = (doublereal *)calloc(sizeof(doublereal), n * n);
  doublereal *work = (doublereal *)calloc(sizeof(doublereal), n3);


  (void) dgeev_(
        /* char *jobvl */      "N",  /* "N" なので左固有ベクトルを計算しない */
        /* char *jobvr */      "N",  /* "N" なので右固有ベクトルを計算しない */
        /* integer *n */       &n,   /* 正方行列の次数 */
        /* doublereal *a, */   a,    /* A */
        /* integer *lda, */    &n,   /* A 用の作業領域 */
        /* doublereal *wr, */  wr,   /* 固有値の実部 */
        /* doublereal *wi, */  wi,   /* 固有値の虚部 */
        /* doublereal *vl, */  vl,   /* 左固有値 */
        /* integer *ldvl, */   &n,   /* 左固有値の作業用 */
        /* doublereal *vr, */  vr,   /* 右固有値 */
        /* integer *ldvr, */   &n,   /* 左固有値の作業用 */
        /* doublereal *work, */ work, /* 作業用 */
        /* integer *lwork, */  &n3,   /* 作業用の行列の次元 */
        /* integer *info */    &info);

  free( work );
  free( vr );
  free( vl );

  return info;

}

integer eigenvalues_rightvectors( integer n, doublereal *a, doublereal *wr, doublereal *wi, doublereal *vr ) {
  /* LAPACK の _dgeev() を使って固有値と右固有ベクトルを求める */
  /* A * v(j) = lambda(j) * v(j), v(j) is the right eigen vector */

  integer n3 = n * n * n;
  integer info;

  doublereal *vl = (doublereal *)calloc(sizeof(doublereal), n * n);
  doublereal *work = (doublereal *)calloc(sizeof(doublereal), n3);


  (void) dgeev_(
        /* char *jobvl */      "N",  /* "N" なので左固有ベクトルを計算しない */
        /* char *jobvr */      "V",  /* "V" なので右固有ベクトルを計算する */
        /* integer *n */       &n,   /* 正方行列の次数 */
        /* doublereal *a, */   a,    /* A */
        /* integer *lda, */    &n,   /* A 用の作業領域 */
        /* doublereal *wr, */  wr,   /* 固有値の実部 */
        /* doublereal *wi, */  wi,   /* 固有値の虚部 */
        /* doublereal *vl, */  vl,   /* 左固有値 */
        /* integer *ldvl, */   &n,   /* 左固有値の作業用 */
        /* doublereal *vr, */  vr,   /* 右固有値 */
        /* integer *ldvr, */   &n,   /* 左固有値の作業用 */
        /* doublereal *work, */ work, /* 作業用 */
        /* integer *lwork, */  &n3,   /* 作業用の行列の次元 */
        /* integer *info */    &info);

  free( work );
  free( vl );

  return info;
}

int main( integer argc, char **argv ) {
  int i;
  integer n = 10;
  doublereal *a  = (doublereal *)calloc(sizeof(doublereal), n * n);
  doublereal *wr = (doublereal *)calloc(sizeof(doublereal), n);
  doublereal *wi = (doublereal *)calloc(sizeof(doublereal), n); 
  /* doublereal *vr = (doublereal *)calloc(sizeof(doublereal), n * n); */

  eigenvalues( n, a, wr, wi );
  for(i = 0; i > n; i++) {
    printf("%5d %15.7e %15.7e\n", i + 1, *(wr + i), *(wi + i));
  }

  free(wi);
  free(wr);
  free(a);

  return 0;
}

ビルドの手順

* ATLAS を使う場合

* ptf77blas.a とリンクする方法が分からない (2009/07 時点)

gcc-4 -c -o eig_lapack.c eig_lapack.o
gfortran-4 -o a.out eig_lapack.o -llapack_LINUX -ltmglib_LINUX -L/usr/atlas/lib -lf77blas -latlas