2020年3月 – うだうだWeblog2

2020-03-252020-03-25

LAPACK 3.9.0をGCCでインストール

　気が付いたらLAPACK 3.6.0のインストール方法を書いてから大分経ってたので，Ubuntu環境にTAR ballからインスト―する方法をメモ書きしておく。前提として

GCC (gcc, gfortran)が使える
make (GNU make)が使える
/usr/local/libにライブラリファイル(*.a)を，/usr/local/includeにインクルードファイル(*.h)をインストールできる権限(sudo)がある

を仮定している。 $ 太字がコマンドである。ついでに拙著「LAPACK/BLAS入門」のサンプルプログラム実行方法も書いておいた。

1. LAPACK 3.9.0のTAR ball(v3.9.0.tar.gz)をダウンロード

$ wget https://github.com/Reference-LAPACK/lapack/archive/v3.9.0.tar.gz
–2020-03-25 18:38:15– https://github.com/Reference-LAPACK/lapack/archive/v3.9.0.tar.gz
Resolving github.com… 52.192.72.89
Connecting to github.com|52.192.72.89|:443… connected.
HTTP request sent, awaiting response… 302 Found
Location: https://codeload.github.com/Reference-LAPACK/lapack/tar.gz/v3.9.0 [following]
–2020-03-25 18:38:16– https://codeload.github.com/Reference-LAPACK/lapack/tar.gz/v3.9.0
Resolving codeload.github.com… 52.193.111.178
Connecting to codeload.github.com|52.193.111.178|:443… connected.
HTTP request sent, awaiting response… 200 OK
Length: unspecified [application/x-gzip]
Saving to: “v3.9.0.tar.gz”
[ <=> ] 7,534,567 4.84M/s in 1.5s
2020-03-25 18:38:18 (4.84 MB/s) – “v3.9.0.tar.gz” saved [7534567]
$ ls -l v3.9.0.tar.gz
-rw-rw-r– 1 tkouya tkouya 7534567 Mar 25 18:38 v3.9.0.tar.gz
$ tar zxvf v3.9.0.tar.gz
lapack-3.9.0/
lapack-3.9.0/.appveyor.yml
(略)
lapack-3.9.0/meson.build
lapack-3.9.0/meson_options.txt
$

2. “lapack-3.9.0″ディレクトリに移動してMakefile用のincludeファイル(make.inc)をgfortran用に設定

$ cd lapack-3.9.0
$ ls
BLAS CTestCustom.cmake.in lapack.pc.in meson.build
CBLAS DOCS lapack_testing.py meson_options.txt
CMAKE INSTALL LICENSE README.md
CMakeLists.txt lapack_build.cmake Makefile SRC
CTestConfig.cmake LAPACKE make.inc.example TESTING
$ ls -l INSTALL
total 228
-rw-rw-r– 1 tkouya tkouya 483 Nov 21 16:57 CMakeLists.txt
-rw-rw-r– 1 tkouya tkouya 5355 Nov 21 16:57 dlamch.f
（略)
-rw-rw-r– 1 tkouya tkouya 2140 Nov 21 16:57 slamchtst.f
-rw-rw-r– 1 tkouya tkouya 21748 Nov 21 16:57 tstiee.f
$ cp INSTALL/make.inc.gfortran ./make.inc

3. Reference BLASを生成

$ cd BLAS
$ ls
blas.pc.in CMakeLists.txt Makefile SRC TESTING
$ make
make -C SRC
(略)

4. CBLASを生成

$ cd ../
$ cd CBLAS
$ ls
cblas.pc.in CMakeLists.txt include README testing
cmake examples Makefile src
$ make
cp include/cblas_mangling_with_flags.h.in include/cblas_mangling.h
make -C src

5. librefblas.a, libcblas.aが生成されていることを確認

$ cd ../
$ ls -l
-rw-rw-r– 1 tkouya tkouya 384424 Mar 25 18:44 libcblas.a <– cblas
-rw-rw-r– 1 tkouya tkouya 618312 Mar 25 18:43 librefblas.a <– reference BLAS

6. LAPACKを生成

$ make
make -C INSTALL run
make[1]: Entering directory `/home/tkouya/pool/lapack/lapack-3.9.0/INSTALL’
gfortran -O2 -frecursive -c -o lsame.o lsame.f
（略)
–> LAPACK TESTING SUMMARY <–
Processing LAPACK Testing output found in the TESTING directory
SUMMARY nb test run numerical error other error
================ =========== ================= ================
REAL 1308195 0 (0.000%) 0 (0.000%)
DOUBLE PRECISION 1309007 0 (0.000%) 0 (0.000%)
COMPLEX 760590 1 (0.000%) 0 (0.000%)
COMPLEX16 768086 1 (0.000%) 1 (0.000%)

–> ALL PRECISIONS 4145878 2 (0.000%) 1 (0.000%)

7. LAPACKEを生成

$ cd LAPACKE
$ make
(略)
make[1]: Leaving directory `/home/tkouya/pool/lapack/lapack-3.9.0/LAPACKE/utils’

8. ファイル確認

$ cd ../
$ ls -l
-rw-rw-r– 1 tkouya tkouya 384424 Mar 25 18:44 libcblas.a
-rw-rw-r– 1 tkouya tkouya 12336044 Mar 25 18:51 liblapack.a
-rw-rw-r– 1 tkouya tkouya 6898692 Mar 25 19:11 liblapacke.a
-rw-rw-r– 1 tkouya tkouya 618312 Mar 25 18:43 librefblas.a
-rw-rw-r– 1 tkouya tkouya 630092 Mar 25 18:51 libtmglib.a

9. /usr/local/libにライブラリ(*.a)ファイルを，/usr/local/includeにヘッダファイル(*.h)をコピー

$ sudo cp *.a /usr/local/lib
$ sudo cp CBLAS/include/*.h /usr/local/include
$ sudo cp LAPACKE/include/*.h /usr/local/include
$ ls -l /usr/local/lib
-rw-r–r– 1 root root 384424 Mar 25 19:15 libcblas.a
-rw-r–r– 1 root root 12336044 Mar 25 19:15 liblapack.a
-rw-r–r– 1 root root 6898692 Mar 25 19:15 liblapacke.a
-rw-r–r– 1 root root 618312 Mar 25 19:15 librefblas.a
-rw-r–r– 1 root root 630092 Mar 25 19:15 libtmglib.a
$ ls -l /usr/local/include
-rw-r–r– 1 root root 21243 Mar 25 20:42 cblas_f77.h
-rw-r–r– 1 root root 29918 Mar 25 20:42 cblas.h
-rw-r–r– 1 root root 444 Mar 25 20:42 cblas_mangling.h
-rw-r–r– 1 root root 8051 Mar 25 20:42 cblas_test.h
-rw-r–r– 1 root root 4236 Mar 25 19:15 lapacke_config.h
-rw-r–r– 1 root root 834528 Mar 25 19:15 lapacke.h
-rw-r–r– 1 root root 474 Mar 25 19:15 lapacke_mangling.h
-rw-r–r– 1 root root 33161 Mar 25 19:15 lapacke_utils.h
-rw-r–r– 1 root root 455216 Mar 25 19:15 lapack.h

10. 「LAPACK/BLAS入門」サンプルプログラムをインストール

$ git clone https://github.com/tkouya/lapack_blas_tutorial
$ cd lapack_blas_tutorial

Makefile.unix → 下記のように書き換え

#*************************************************#
# LAPACK/BLAS Tutorial                            #
# Makefile for Linux gcc or icc environment       #
# Last Update: 2016-12-02 (Fri) T.Kouya           #
#*************************************************#
# ------------------------ start of configuration --------------------------
# Intel C compiler
#include lapack_icc.inc ←コメント化

# GNU Compiler Collection
include lapack_gcc.inc ← コメント外して有効化

（略)
# ------------------------ end of configuration --------------------------
all: first matvec_mul linear_eq blas row_column_major eig dgecon openmp pthread ← "gpu"を外す
#all: first matvec_mul linear_eq blas row_column_major gpu eig dgecon openmp pthread spblas_mkl

lapack_gcc.inc → 下記のように書き換え

#*************************************************#
# LAPACK/BLAS Tutorial                            #
# Configuration file for GNU compiler collection  #
# Last Update: 2016-12-02 (Fri) T.Kouya           #
#*************************************************#
CC=gcc
FC=gfortran
CPP=g++

INC = -I/usr/local/include
#LIB = -L/usr/local/lib/gcc  -lgfortran -lm
LIB = -L/usr/local/lib  -lgfortran -lm

CBLAS_INC = $(INC)
CBLAS_LIB = $(LIB) -lcblas -lrefblas -lgfortran

LAPACKE_INC = -I/usr/local/include/lapacke $(CBLAS_INC)
#LAPACKE_LIB = -L/usr/local/lib/gcc -llapacke -llapack $(CBLAS_LIB)
LAPACKE_LIB = -L/usr/local/lib -llapacke -llapack $(CBLAS_LIB)

#IMKL_INC=-I/opt/intel/mkl/include
#IMKL_LIB=-L/opt/intel/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -L/opt/intel/lib/intel64 -lifcore
#IMKL_LIB=-L/opt/intel/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -L/opt/intel/lib/intel64 -lifcore -liomp5

OPENMP = -fopenmp
OPENMP_LIB = -lgomp

Makeしてコンパイル→実行ファイルの一つをお試し

$ make -f Makefile.unix
$./lapack_dgeev ←行列の固有値・固有ベクトルを計算
Dim = 3 ←次元数を入力
Eigenvalues =
0: ( 1.82216, 0)
1: ( 0.103726, 0.366921)
2: ( 0.103726, -0.366921)

(0.619085, 0) (0.634996, 0) (0.634996, 0)
(0.641249, 0) (-0.399111, -0.289681) (-0.399111, 0.289681)
(0.453358, 0) (-0.396179, 0.443416) (-0.396179, -0.443416)
||A * right_ev – lambda * right_ev||_F / ||A||_F = 4.35690729441612009e-16
||A^T * left_ev – lambda * left_ev||_F / ||A||_F = 3.79955740713784858e-16

2020-03-162020-03-16

3/16(月) 袋井・晴

　寒の戻りとかで，２月中旬並みの気温，東京では雪まじりのみぞれが降ったらしい。こちら静岡も天気は良いが手袋・マフラー・ニット帽の完全防備体制でもチト寒いぐらい。まぁすぐに治るようなので，基準木が開花するのを楽しみに待つことにするのである。

　先週土曜日は本来本学の卒業式だったのだが，コロナウイルス騒動の煽りで中止，市内で感染者が出たこともあり，小規模の手交式も中止となったのを憂いた若手教員の方々が「オンライン卒業式」を開催して頂き，ワシも額縁ょ〜として挨拶を執り行った。状況が状況だけに，タイムリーなローカルニュースネタとなり，いろいろな方々に「観たよ〜」というご連絡を頂いた次第。お声がけいただいた方々，ご準備頂いた情報学部の先生方に感謝すると共に，先週木曜日に散髪を済ませておいて良かったと神さん共々安堵したのでありました。詳細は本学blog記事をご覧頂きたい。動画はNHKとかに上がってたらしいので適当に探して下さいまし。

　今日はその代休日なんだけど，静かでいいかと思ってきてみれば，計画停電後で学内LANのSwitchが全然起動してないでやんの。外部はおろか，内部の別の階のマシンにも繋がらない（午後3時ごろから復活しました。休日出勤ありがとうございます）。つ〜ことでこれは自分のギガ消費しながらテザリングで書いてますのでこの辺で。少しはPythonコード書かなきゃなぁ。

2020-03-112020-03-12

コロナウィルス流行曲線をPythonで描いてみる

　佐藤總夫「自然の数理と社会の数理 II 微分方程式で解析する」（日本評論社）という名著があった。ホントはI, II, IIIと３分冊で出るはずが，I, IIが出た後は音沙汰なしという，当時学部生だったワシが残念に思っていたシリーズである。今回のコロナウィルス騒ぎで散々出てくるあの「流行曲線」という奴，本書に詳しく出てたので，折角だから知識と本の虫干しがてら，Python+SciPyで描いてみることにしよう。つーことで，以下の数学的な記述は全部佐藤著の「第７話　伝染病の伝播」に基づくものなので，詳しく知りたい人は図書館にGoである。理工系大学ならあるんじゃないのかな。

　感染症モデルの中でも一番簡単なものを取り上げる。前提として

感染者を[latex]y = y(u)[/latex], 感受性者（これから感染する可能性のある人）を$x(u) = x$とする。
全人口を$n+1$とし，増減はないものとする（「死者は出ない」の意）。つまり常に$x + y = n + 1$が保持される。
感染者の発生率は$xy$に比例する（人々の均等な接触を表現）。
最初の感染者は$y(0) = 1$とする。

とする。この時，常微分方程式の初期値問題としては

\[ \frac{dx}{du} = -xy = -x(n+1 – x), x(0) = n \]

となる。新たにどの程度の感染者が現れるかを表現したもの，つまり感染スピードを表わす関数が$w(u)$である。上記を解いて$x(u)$, $y(u)$を求め，そこから$w(u)$を計算してグラフ化したものが下記の図である。小うるさい人にすれば，人数が単位の$x(u)$, $y(u)$と人/時間が単位の$w(u)$を同じ縦軸に描くとは何事だと言いそうであるが，これらの関数の関係を表したものとご理解頂きたい・・・今気がついたが，右軸も使えば良かったんだな，面倒くさいから，やりたい人向けの課題としておく。

　このグラフから，感覚的には次のことはすぐに分かる。

感染を妨げるものはないので，最終的には全員に感染する。
感染が止まるのは全員に感染するからであって，感染するスピードは山なりに変化する。

　初期感染者は今更ゼロにはできないし，できることといえば，感染者を見つけたら素早く感受性者から引き離して$xy$ を減らすしかない。つまりは初期感染スピード$w(0)$を極力減らすしかない訳で，実際，$w(0) = 10, 20, 30, 40$まで$10$刻みで変えてみると下記のようなグラフを得る。

　このグラフが感染スピードの変化を示しており，当然ながら，初期感染スピードが遅ければ遅いほど感染率極力小さく抑えられることが分かる。現在の厚労省が取っている対策もこれに基づくもので，このグラフはまさしくこのグラフが描いている曲線ということになる。

　つーことで，上のグラフを描くためのPython + SciPyコード例を下記に挙げておく。ちなみに，この常微分方程式は解析解が手計算で求められるので，ピーク時はどこになるとかといった議論も可能である。詳細は佐藤著をご覧頂きたい。

# ode_epidemic.py: 常微分方程式の初期値問題

# numpy
import numpy as np

# integrateパッケージ
import scipy.integrate as spint

# matplotlib
import matplotlib.pyplot as plt

num_div_t = 100
windex = 0

w = [[0.0] * num_div_t for i in range(4)]
#print('w[0] = ', w[0])

# 人口
all_num_people = [11, 21, 31, 41]

for num_people in all_num_people:
    # 初期感染者
    num_infected_people = 1
    # 初期感受性者
    num_noninfected_people = num_people - num_infected_people

    # 陽的形式の右辺
    # y' = func(t, x) = -x(n+1-x)
    def func(t, y):
    #	return [-y[0] * (num_people - y[0])]
        return np.array([-y[0] * (num_people - y[0]), y[1] * (num_people - y[1])])

    # 初期値
    # y(0) = [n, 1]
    #y0 = [num_noninfected_people]
    y0 = [num_noninfected_people, num_infected_people]

    # t = [0, 1]
    t_interval = [0.0, 1.0]
    #print(t_interval)

    # 常微分方程式を解く
    #ret = spint.solve_ivp(func, t_interval, y0)
    ret = spint.solve_ivp(func, t_interval, y0, t_eval = np.linspace(t_interval[0], t_interval[1], num_div_t))

    # 感染者の増加率を計算
    w[windex][0] = num_noninfected_people # 初期感染者発生率
    for index in range(1, ret.t.size):
        w[windex][index] = (ret.y[1, index] - ret.y[1, index - 1]) / (ret.t[index] - ret.t[index - 1])
    windex += 1

    # 結果を表示
    #print(ret)

    # yを表示
    #print(ret.y)

# wを表示
#print('w[0] = ', w[0])

# t-yグラフを描画
# グラフ初期化
figure, axis = plt.subplots()

# 値をセット
#axis.plot(ret.t, ret.y[0, :], label='x(u)')
#axis.plot(ret.t, ret.y[1, :], label='y(u)')
axis.plot(ret.t, w[0], label='w(0) = '+ str(all_num_people[0] - 1))
axis.plot(ret.t, w[1], label='w(0) = '+ str(all_num_people[1] - 1))
axis.plot(ret.t, w[2], label='w(0) = '+ str(all_num_people[2] - 1))
axis.plot(ret.t, w[3], label='w(0) = '+ str(all_num_people[3] - 1))

# x軸，y軸，グラフタイトルをセット
axis.set(xlabel = 't', ylabel = 'w', title = 'Epidemic model')

# グリッドを描画
axis.grid()

# 凡例
axis.legend()

# グラフ保存ファイル名
figure.savefig('ode_epidemic_w.png')

# グラフを画面に描画
plt.show()

2020-03-082020-03-08

3/8(日) 駿府・雨

　山松ゆうきちの漫画のような「さかざんざんざん」という土砂降りの日。冬らしくないが，春を感じさせるような暖かいものでもない氷雨。今年はソメイヨシノの開花が次週にも迫っているようだが，そう簡単にいくかどうか。次週末は駿府城公園を探索してみよう。標準木，どうなっているかなぁ。

　静岡市内の警察署に出頭してきた。免許更新なのだが，ワシは優良ドライバーなので，圏内のどこの警察署でも30分のビデオを見るだけで更新できるのである。うっかりして床屋に行きそびれてしまい，髪も髭も伸び放題，ボサボサ状態の写真で写ることになった。ま，５年間の辛抱である。

　本学紀要原稿を書き上げ，それを流用してPython数値計算入門（仮）の多倍長精度計算の章を全面的に書き直した。今のところPythonにはQDのようなマルチコンポーネント型多倍長精度演算の高速なライブラリがないようなので，紀要原稿ではワシのCライブラリベースのモノを使ってのベンチマークテストを行い，Python数値計算の方はmpmathとgmpy2の比較を詳細に行った。ちと手間であったが，還暦まではせめて紀要原稿の一本ぐらいはコンスタントに書いていきたいものである。
　コロナウイルスが蔓延ってあらゆる出張やイベント予定がパーになっている今月をチャンスとし，Python本の方は偏微分方程式の章の書き足しと，補間計算の解説全面改稿，疎行列とソルバーの記述追記を目指したい。これができれば，あとはグラフィックスの方をどこまで追記できるかだなぁ。頑張りましょう。

　神さんとコロナ騒ぎで閑散とした街中を散策するも，ラーメン屋には客が入っているし，人通りは少ないが，人々の外出意欲が全くなくなったということではなさそう。マスクもあるところにはあるモンで，まだ我が家には備蓄はあるが，少し買い足しておいた。どこに売っていたかは書いてあげないのである。ワシらは老舗の洋食屋さんで昭和チックなランチセットを食して満足した。

　そーいや，ボチボチ額縁ょ〜就任一年になるんだなぁ。と言っても慣らし運転もいいところで，未だにあれこれ迷うこと多かれど，色々な情報が集まってくると学習するところも多かった。今後の大仕事は教員評価だが，大筋，自分が見知っているところと違う点は少ないことも分かった。まぁ査読論文・査読付き国際研究集会・著作等々，第三者の目が通ってパブリッシュされたものがどんだけあるかで仕事ぶりは分かる訳で，どんだけ言い訳しようともそれが全て。「論文書いて」と言えば済む，簡単なお仕事です額縁ょ〜は。任期はあと一年なので，せいぜい小うるさく嫌がられながら頑張りましょう。

　次週はPython仕事しながら今月末に締め切りのある国際研究集会用の仕事せにゃ。コロナ騒ぎのグローバル化でどこの集会の開催も危ういけど，人にハッパかけるなら自分もある程度はやらんとねぇ。声のでかい役立たずの評論家モドキにはなるまいぞなるまいぞ。

　風呂入って寝ます。