このサイトについて・著者連絡先
2020/03/08
とある必要性のため疎行列計算ライブラリを調べていて、SuperLU というのを2014年の暮れ頃に見つけました。
以来、たまーに気が向いたときに動向を覗いています。
本家ポータルサイトはここ。
GitHub には
逐次版、
マルチスレッド (OepnMP) 版、
分散メモリー版
（以上最新のリンクに更新）がありますが、以下主に逐次版について述べます。
逐次版5.2 および 5.3.0、MS-VC++2013 (Express) で確認。意外と簡単です。わざわざ Makefile を書き換えたり統合開発環境に取り込んだりせずともこれで十分です。ただし、以下ではわざと、警告メッセージをできるだけ表示させる設定としています（ソースが正しく記述されていれば実害がほぼない C4210（非標準の拡張機能）と C4996（非セキュア版 CRT 関数）のみを無効化）
superlu_5.2.0\CBLAS に移動。
cl /c /O2 /MD /W4 /wd4210 *.c
lib /out:cblas.lib *.obj
superlu_5.2.0\SRC に移動。
cl /c /O2 /MD /W4 /wd4210 /wd4996 *.c
lib /out:superlu.lib *.obj
これを実行するとおびただしい数の警告が出されるとおり、ソースコード全般にちょっと詰めが甘いところが見受けられます。32bit, 64bit ともに、一度は警告メッセージにつぶさに目を通し、実際に使う前に、SuperLU のもつ潜在的な問題点を正しく認識・掌握しておくことをお勧めします。
成果物は superlu.lib と cblas.lib です。これらを使用するプログラムのリンク時に追加します。
細々とした修正が必要ですが、ビルドできることは判りました。執筆中。
小規模な計算ではかえって遅くなります。
言うまでもなくMinimalist GNU for Windowsでのビルドは可能で、もちろんちゃんと動きます。
なんでこんな古いバージョンのGCCを使うのか？と言われそうですが、ちょっと複雑な事情がありまして……。あくまで事例紹介です。普通の用途なら新しいGCCや、Clangを使うのが良いでしょう。
make.inc があります。これを、以下の定義が有効になるように編集します（重複する既存定義はコメントアウト）。
SuperLUroot = ..
SUPERLULIB = $(SuperLUroot)/lib/libsuperlu_5.0.a
BLASLIB = $(SuperLUroot)/lib/libblas.a
BLASDEF = -DUSE_VENDOR_BLAS の既存定義をコメントアウト無効化する（Windowsにはベンダー提供BLASライブラリは当然無く、SuperLU同梱のCBLASを使うため）。
SuperLUroot の定義は「SRCディレクトリから見て一つ上」の意味、すなわち、書庫を展開した最上位ディレクトリを指します。libsuperlu_5.0.a は、バージョン4.3なら libsuperlu_4.3.a 等とすれば良いです。
cd する。
make blaslib を実行してCBLASをビルド。makeは、このようにタイプするとMSYSの make.exe が使われます。当然、MinGW の mingw32-make.exe を使っても同じ結果になります。
make を実行してSuperLU本体をビルド。
lib ディレクトリに、出来上がった libblas.a と libsuperlu_5.0.a があります。後はこれをunix風にgccのlibにコピーしても良いし、その都度指定しても良い。
Makefileの中身をいじる（unix由来コマンドを対応するDOSコマンドに変更し、マルチステートメントをシングルステートメントに分解する等）と、MSYSなしでも何とかビルドすることもできるのですが、まあ面倒です。素直にMSYSを使う方が作業効率は良いです。
使い方は以下のようなやりかたなどとします。
set LIBRARY_PATH=[libsuperlu_5.0.a と libblas.a を置いた場所];%LIBRARY_PATH%
g77 -O2 foo.for -lsuperlu_5.0 -lblas -lsuperlu_5.0
バージョン5.0をgcc-3.4.5でビルドすると、なぜか、このように -lblas を挟んで重複指定をしないと「参照が解決できない」と言われてしまいます。バージョン4.3では１回指定で問題ありません。
SuperLU_MT 3.0もgcc-3.4.5でビルドできることはできますが、gcc-3ではpthread限定かつオーバーヘッドがかなり大きく、速度メリットを得にくい模様です。
SuperLUを見つけたのは、Pythonの疎行列計算ライブラリ Pysparse がこれに依存していると書いてあったからです。SuperLUのライセンスは3条項BSD（修正BSD）であり、今回の調査の目的に適うものでした。
ライセンスについて注意喚起
本体部分は基本的には修正BSDだが、同梱されているいくつかのソースコードに、ライセンスが異なる他パッケージからの流用物があり、取り扱いが不適切であるとして Debian では議論になっている模様です。https://lists.debian.org/debian-science/2016/05/msg00115.html
ちょっと嫌な予感がしていたのですが、やはりといった感じです。当面、動向を注視します。
ユーザー責任でよく見て使用する必要がありそうです。
Pysparseは他にもSuiteSparse (http://faculty.cse.tamu.edu/davis/suitesparse.html) に含まれる UMFPACK も採用していて、使用者が任意に選択（併用）できるようです。こちらは GNU GPL2 or later とのことで、LGPLですらないため残念ながら目的にそぐわず。
SuperLU, UMFPACKとも、「京」などにも移植されているそうです。
サンプルプログラム superlu.c において、求解後のB（解が上書きされている）を表示させてみる：
dPrint_Dense_Matrix("B", &B);
Dense matrix B: Stype 6, Dtype 1, Mtype 0 nrow 5, ncol 1, lda 5 nzval: -0.031250 0.065476 0.013393 0.062500 0.032738
一方、サンプルプログラムと同じ方程式をノーマルのNumPy（内部的にはReference BLASとLAPACK）でやると
>>> import numpy as np
>>> s=19
>>> u=21
>>> p=16
>>> e=5
>>> r=18
>>> l=12
>>> A=np.array([
... [s, 0, u, u, 0],
... [l, u, 0, 0, 0],
... [0, l, p, 0, 0],
... [0, 0, 0, e, u],
... [l, l, 0, 0, r],
... ], dtype=np.float64)
>>> B=np.array([1, 1, 1, 1, 1], dtype=np.float64)
>>> A
array([[ 19., 0., 21., 21., 0.],
[ 12., 21., 0., 0., 0.],
[ 0., 12., 16., 0., 0.],
[ 0., 0., 0., 5., 21.],
[ 12., 12., 0., 0., 18.]])
>>> B
array([ 1., 1., 1., 1., 1.])
>>> np.linalg.solve(A,B)
array([-0.03125 , 0.06547619, 0.01339286, 0.0625 , 0.0327381 ])
単に printf() の書式指定子の違いが出ているだけであろう。
[2020/02/29] 付属の CBLAS はコンパイラが未初期化変数警告を出してくるなど、ちょっと怪しいところがあるようです。
そこで、Reference BLAS の最新版（これを書いている 2020/02/29 現在 3.8.0）を f2c 経由でビルドする方法 を調べてみました。付属 CBLAS と差し替えたところ、dgssv() だけしか試していないがとりあえずはうまく動くようです。
[2023/05/22] 2023年4月に SuperLU 6.0.0 が出ているようです。Changelog に "Add 64-bit indexing support, with mixed int and int64_t types." という記載がありますので、下記はもう古い情報かもしれませんが詳細未調査です。おそらくC言語からの利用に限定で、FORTRAN77 ではその言語仕様に size_t の概念が存在しないためどのみちダメでしょう。
[2020/03/08] ポインタの差を単なる int に代入している個所がいくつかあり、LLP64 (Windows における 64bit ビルド) では MSVC の警告レベル /W3 以上で次の警告が出ます。
warning C4244: '=' : '__int64' から 'int' への変換です。データが失われる可能性があります。
これは最悪の場合ポインタ破損につながる問題であり、実用上も問題が生じる可能性があるかもしれません。32bit および LP64 なら問題ないと思われるが、さて？ （なお、C4244 は double から float への変換、int から float への変換などでも多数出てくるが、ここでは本当に不味いポインタ破損のみについて議論）
LLP64 データ型モデルでは、たぶん、まずい箇所を intptr_t に変更する手当てが必要なのだと思います。
場所によっては（ extra 、top1 、top2 ）意味的に ptrdiff_t のほうがより正しいかもしれないが、MSVC でターゲットアーキテクチャを amd64 とすると intptr_t と ptrdiff_t はいずれも内部的に __int64 で typedef されることになり実質可換なので、ひとまずこれで。要 C99/C++11。
逐次版 version 5.2.1 からの差分は以下の通り（逐次版 ver.5.3.0 でもこの LLP64 問題は放置されたままで、同じ修正が必要なようです）。
なお、この修正を施したとしても、多くの関数で要素数やそれに準じるものを int で宣言しているので、LLP64 では結局、32bit のキャパシティーを超えるものは扱えないことに注意。FORTRAN77 インターフェイスを持っている以上これは仕方がないことです（FORTRAN77 には size_t の概念が存在しない）。
memory.c
132c132
< user_bcopy(char *src, char *dest, int bytes)
---
> user_bcopy(char *src, char *dest, intptr_t bytes)
cmemory.c
35c35
< extern void user_bcopy (char *, char *, int);
---
> extern void user_bcopy (char *, char *, intptr_t);
342c342,343
< int isize, dsize, extra;
---
> int isize, dsize;
> intptr_t extra;
502c503,504
< int new_len, tries, lword, extra, bytes_to_copy;
---
> int new_len, tries, lword;
> intptr_t extra, bytes_to_copy;
dmemory.c
35c35
< extern void user_bcopy (char *, char *, int);
---
> extern void user_bcopy (char *, char *, intptr_t);
342c342,343
< int isize, dsize, extra;
---
> int isize, dsize;
> intptr_t extra;
502c503,504
< int new_len, tries, lword, extra, bytes_to_copy;
---
> int new_len, tries, lword;
> intptr_t extra, bytes_to_copy;
smemory.c
35c35
< extern void user_bcopy (char *, char *, int);
---
> extern void user_bcopy (char *, char *, intptr_t);
342c342,343
< int isize, dsize, extra;
---
> int isize, dsize;
> intptr_t extra;
502c503,504
< int new_len, tries, lword, extra, bytes_to_copy;
---
> int new_len, tries, lword;
> intptr_t extra, bytes_to_copy;
zmemory.c
35c35
< extern void user_bcopy (char *, char *, int);
---
> extern void user_bcopy (char *, char *, intptr_t);
342c342,343
< int isize, dsize, extra;
---
> int isize, dsize;
> intptr_t extra;
502c503,504
< int new_len, tries, lword, extra, bytes_to_copy;
---
> int new_len, tries, lword;
> intptr_t extra, bytes_to_copy;
slu_util.h
310c310
< int size; /* length of the memory that has been used */
---
> intptr_t size; /* length of the memory that has been used */
315,318c315,318
< int size;
< int used;
< int top1; /* grow upward, relative to &array[0] */
< int top2; /* grow downward */
---
> intptr_t size;
> intptr_t used;
> intptr_t top1; /* grow upward, relative to &array[0] */
> intptr_t top2; /* grow downward */
この文書が誰か・何かの役に立てば幸いです。さらに願わくば、この文書の公開（並びに Lis / Eigen / MUMPS の各事例）が、日本における数値計算ライブラリ利用裾野を広げ、日本全体の底上げにつながることを、静かに期待します。
トップページ /
計算機関係の雑多なこと /
SuperLUをWin32上でビルドする
このサイトについて・著者連絡先