ノート/EM?
訪問者数 5966      最終更新 2012-07-23 (月) 16:13:55

LAPACK (Linear Algebra Package)とは (2012-07-13)

ウィキペディア LAPACK 曰く

LAPACK (Linear Algebra PACKage) は線型計算のための数値解析ソフトウェアライブラリで、
線型方程式や線型最小二乗問題、固有値問題、特異値問題等を数値的に解くために利用される。
本ライブラリは複素数または実数を成分とする行列を扱うことが可能であり、LU分解やコレスキー分解、
QR分解、シュア分解(英語版)等の行列分解(英語版)を行うためのサブルーチンを含む。
サブルーチンは単精度版と倍精度版が提供される。1992年のLAPACKの初版はFORTRAN 77で実装されていたが、
現在はFortran 90が用いられている。LAPACK 3.4.0からはC言語インターフェースである LAPACKEが統合され、
C言語やC++からの利用が容易になった。

LAPACKはLINPACKおよびEISPACKの後継と見做されている。

というわけで、この辺がよさそう。

元のFortran版のLAPACKは、netlib.orgのlapackのページ にある。 最新はLAPACK 3.4.1 (2012/4/20付)。

いろいろな「派生」(と呼んでいいかどうか疑問だが)があるようで、たとえばここのFAQにいろいろと書いてある。

非常に短く言うと、

というわけで、LAPACKの本体パッケージ+LAPACKE あたりが本命だろうか。

別の問題として、LAPACKを支えるBLASをどうするか、がある。BLASはベクトル・行列の単純な操作をパッケージ化したものだが、この部分をチューニングすることで高速化を図る。Lapackパッケージに付属しているreference implementationは機種依存のチューニングがされていないため、十分に速くはないという。これに対して、自動的に機種依存のチューニングを行うATLAS(上でも参照、かなり汎用的)や、Intel x86系を特に対象としたGoto2、各ハードメーカーが作っている様々なバージョン、などがあるようだ。

そこで、全体としては

ATLAS + LAPACKを考える。
ATLASのパッケージから参照されるLAPACKを、同時にコンパイルすることによって、LAPACKから(LAPACKに付いているBLAS reference implementationではなくて)ATLASを使うようにする
このモードだと、C言語インターフェースであるLAPACKEは含まれないので、別途LAPACKパッケージ中のLAPACKE部分をコンパイルしライブラリ化する。(但し、手抜きをしたのでLAPACKE部分だけはダイナミックリンクライブラリではなくデフォルトのスタティックリンクになった。手間を掛ければできそうだが。)

ATLAS + LAPACKのインストール

yumで見ると、手元のマシンではLAPACKの最新rpmパッケージが3.2.1なので、欲しい3.4以上でない。ソースからインストールすることにする。

Lapackの単独インストール(2012-07-18) 〜 意味なかった

というわけで、LAPACKのベースインストールは完了。

ATLASをインストールする時に、一緒にLAPACKをインストール、あとからLAPACKEをインストール(2012-7-20)〜結局これが正しかった

ここまででatlasとlapackのライブラリやincludeファイルができる。

ls /usr/local/atlas/lib
libatlas.a  libf77blas.a  libptf77blas.a  libtatlas.so libcblas.a
liblapack.a   libptcblas.a  libsatlas.so
ls /usr/local/atlas/include
atlas    clapack.h  cblas.h  

LAPACKEは別途makeする必要がある。

使ってみる

次は、自前のプログラムを書いてテストする。といっても良く分からないので、そこらにあったサンプルプログラムをパクることに。

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <assert.h>
#include <time.h>
#include "lapacke.h"

double A[36] =
{
4, 3, 0, 4, 0, 3,
2, 3, 3, 2, 2, 0,
0, 3, 2, 2, 0, 3,
3, 3, 2, 3, 0, 0,
1, 2, 3, 4, 0, 2,
2, 4, 2, 4, 3, 2
};

double Arow[36] =
{
     4,     2,     0,     3,     1,     2,
     3,     3,     3,     3,     2,     4,
     0,     3,     2,     2,     3,     2,
     4,     2,     2,     3,     4,     4,
     0,     2,     0,     0,     0,     3,
     3,     0,     3,     0,     2,     2,
};

double B[24] =
{
2, 2, 3, 2, 4, 0,
3, 1, 3, 2, 1, 0,
2, 4, 3, 4, 2, 0,
1, 4, 2, 4, 0, 4
};

double Brow[24] =
{
     2,     3,     2,     1,
     2,     1,     4,     4,
     3,     3,     3,     2,
     2,     2,     4,     4,
     4,     1,     2,     0,
     0 ,    0,     0,     4,
};

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, char major, int m, int n, double* a, int lda ) {
        int i, j;
        printf( "\n %s (%c major)\n", desc,major );
        if ( major == 'C')
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", a[i+j*lda] );
                printf( "\n" );
        }

        else {
        printf("m=%d,n=%d\n",m,n);
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", a[j+i*lda] );
                printf( "\n" );
        }
        }

}

void eq_solve(double *a, double *b, int n, int nrhs, char major)
{
  int info, i;
  //char trans = 'n';
  lapack_int lda, ldb, *ipiv;
  lda = n;
  ipiv = malloc(sizeof(int)*n);
        /* Print A */
  print_matrix( "A", major, n, n, a, lda );

  if (major == 'C'){
  ldb = n;
        /* Print B */
  print_matrix( "B", major, n, nrhs, b, ldb );
  //EQUIVALENT STATEMENT
  info = LAPACKE_dgesv(LAPACK_COL_MAJOR, (lapack_int) n, (lapack_int) nrhs, a, lda, ipiv, b, ldb);
  //dgesv_( &n, &nrhs, a, &lda, ipiv, b, &ldb, &info );
  }
  else {
  ldb = nrhs;
        /* Print B */
  print_matrix( "B", major, n, nrhs, b, ldb );
  info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, (lapack_int) n, (lapack_int) nrhs, a, lda, ipiv, b, ldb);
  }

  if (info != 0)
  {
    fprintf(stderr, "dgesv: info = %d\n", info);
  }
  assert(info==0);
  printf("dgesv passed\n");
  //debug_print("b using lapacke", b,n,nrhs,0);
  //print answer here
  print_matrix( "Answer", major, n, nrhs, b, ldb );

  free(ipiv);
}

int main(void)
{
  printf("==============================\n");
  printf("\tSOLVING IN COL MAJOR\n");
  printf("==============================\n");
  eq_solve(A,B,6,4,'C');
  printf("\n==============================\n");
  printf("\tSOLVING IN ROW MAJOR\n");
  printf("==============================\n");
  eq_solve(Arow,Brow,6,4,'R');

  return 0;
}

コンパイルのためのMakefileを次のようにした。

CC      = gcc
LIBDIR = -L/usr/local/lapack
INCDIR = -I/usr/local/lapack/lapacke/include
LIBS   = -llapacke -llapack -lrefblas -lgfortran -lm
CFLAGS = -Wall
LDFLAGS =

OBJ    = $(SRC:.c=.o)
SRC    = invtest3.c
TARGET = invtest3

all:    $(TARGET) $(OBJS)

$(TARGET): $(OBJ)
        $(CC) $(LDFLAGS) -o $@ $(OBJ) $(INCDIR) $(LIBDIR) $(LIBS)

.c.o:
        $(CC) -c $< $(CFLAGS) $(INCDIR) $(LIBDIR) $(LIB)

.PHONY: clean
clean:
        rm -f *~ $(TARGET) $(OBJ)

これで計算した結果は

==============================
        SOLVING IN COL MAJOR
==============================

 A (C major)
   4.00   2.00   0.00   3.00   1.00   2.00
   3.00   3.00   3.00   3.00   2.00   4.00
   0.00   3.00   2.00   2.00   3.00   2.00
   4.00   2.00   2.00   3.00   4.00   4.00
   0.00   2.00   0.00   0.00   0.00   3.00
   3.00   0.00   3.00   0.00   2.00   2.00

 B (C major)
   2.00   3.00   2.00   1.00
   2.00   1.00   4.00   4.00
   3.00   3.00   3.00   2.00
   2.00   2.00   4.00   4.00
   4.00   1.00   2.00   0.00
   0.00   0.00   0.00   4.00
dgesv passed

 Answer (C major)
   0.36   0.81  -0.62   0.08
   1.82   2.00  -0.41  -0.53
  -0.69  -0.65  -0.06   0.92
  -1.24  -1.00   1.12   0.29
   0.38   0.77   0.08   0.15
   0.12  -1.00   0.94   0.35

==============================
        SOLVING IN ROW MAJOR
==============================

 A (R major)
m=6,n=6
   4.00   2.00   0.00   3.00   1.00   2.00
   3.00   3.00   3.00   3.00   2.00   4.00
   0.00   3.00   2.00   2.00   3.00   2.00
   4.00   2.00   2.00   3.00   4.00   4.00
   0.00   2.00   0.00   0.00   0.00   3.00
   3.00   0.00   3.00   0.00   2.00   2.00

 B (R major)
m=6,n=4
   2.00   3.00   2.00   1.00
   2.00   1.00   4.00   4.00
   3.00   3.00   3.00   2.00
   2.00   2.00   4.00   4.00
   4.00   1.00   2.00   0.00
   0.00   0.00   0.00   4.00
dgesv passed

 Answer (R major)
m=6,n=4
   0.36   0.81  -0.62   0.08
   1.82   2.00  -0.41  -0.53
  -0.69  -0.65  -0.06   0.92
  -1.24  -1.00   1.12   0.29
   0.38   0.77   0.08   0.15
   0.12  -1.00   0.94   0.35

の通り。

別の人のページ > BLAS, LAPACK チュートリアル パート1  (簡単な使い方とプログラミング)

もう1つ別の人のページ > Linux で GotoBLAS2 バージョン 1.13 と CBLAS と LAPACK バージョン 3.4.1 のダウンロードとビルドとインストール


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-07-23 (月) 16:13:55 (1854d)