数値計算ライブラリの使い方
疎行列向けライブラリを中心に学習
KOBE HPC スプリングスクール（中級）2022 【Zoom開催】
講師：今村俊幸 (理化学研究所計算科学研究センター)

県立大Rokkoサーバでの実行時の注意事項
 環境変数LANGを適切に設定（抜けていることがあるので）% export LANG=C Intelコンパイラ等を利用しますので、ログインノード上で以下を実行して、モジュールを有効にしてください。

% module load intel intelmpi
petsc-intelmpi slepc-intelmpi gnuplot[guest31@rokko1 ~]$ module list
Currently Loaded Modulefiles:
1) intel/19.1.3
4) slepc-intelmpi/3.16.2
2) intelmpi/2020.4.304 5) gnuplot/5.2.7
3) petsc-intelmpi/3.16.2
7-9, March 2022.
Kobe HPC Spring School2ソースファイルとヒントは、~guest31/share/
SLEPc/  SLEPc sample code, a job-script, and Makefile Please ‘cd’ your work directory, then ‘cp -R’ thedirectory (‘%’ is a prompt, and do not type). You willsee
% cp -R ~guest31/share ./third_day% cd third_day
% ls
SLEPc

 他の例題ファイルも多数用意されてます(tests以下を変えて探してください)。

$PETSC_DIR/arch-linux-c-debug/tests/ksp/ksp/tutorials$SLEPC_DIR/arch-linux-c-debug/tests/eps/tutorials
7-9, March 2022. Kobe HPC Spring School3本日の講義 固有値（線形代数）と応用問題

振動問題
 ネットワーク定常問題
 簡単な固有値計算アルゴリズム紹介
 密行列
 疎行列
 一般的な数値計算ライブラリの紹介
 固有値計算ソフトウェアPETSc+SLEPcを使った演習 Rokkoサーバでの実演習







（連立一次方程式）
固有値計算

応用問題への取り組み

7-9, March 2022. Kobe HPC Spring School

4

はじめに

7-9, March 2022. Kobe HPC Spring School5表記法について ベクトル：い）。

小文字で表記（’→’記号はつけない、太文字にもしな
m行(縦方向) 行列：

大文字イタリックで表記する。
m行(縦方向), n列(横方向)の様に, ベクトルを並べた記法使用も 数体：

特に指示がなければ実数であり、複素数は以下の記法を使用する。ただし、は添字としても使用するので、文脈で判断。

7-9, March 2022. Kobe HPC Spring School6

表記法について ノルム：

縦二本線表記（原則２ノルム）

 行列のノルム：

 行列式：

上記のノルムから定義、縦二本線表記

縦１本線表記 条件数：

7-9, March 2022. Kobe HPC Spring School7

表記法について 連立一次方程式

記号を適時変えることがあるが、基本は上記記号を用いることが多い 残差ベクトルをとして定義する。
 （未知の真解との近さを評価できないため）、残差ベクトルのノルムで近似解を評価する。
 固有値固有ベクトル

固有値をギリシャ記号のラムダ( ) もしくはシータ( ) 固有ベクトルはもしくは を用いて表現するが、適時記号を変えるので注意。
 固有値と対応する固有ベクトルを「固有対」と呼ぶ
 固有値の記号に対応した大文字の行列は、固有値を対角に並べてできる対角行列とする
 また、固有ベクトルを並べてできるベクトルも同様に同じ文字の大文字で表記することがある。


7-9, March 2022. Kobe HPC Spring School8固有値(線形代数)と応用問題

主に線形代数の復習から

7-9, March 2022. Kobe HPC Spring School9固有値とは
 工学における現象解析、システム解析に利用する

構造物の耐震振動解析

ネットワークなどの定常状態解析

 計算方法：線形代数の理論上は次の特性多項式を解けばよい



コンピュータ上で行列式の計算は難しい。多項式の求解は？ニュートン法でできるか？

7-9, March 2022. Kobe HPC Spring School10固有値計算の例
 柱の変形

境界条件を考慮して解は7-9, March 2022. Kobe HPC Spring School11多自由度系の振動
 支配方程式系



変数分離法で解ける



行列・ベクトルで整理するとの形をしている7-9, March 2022. Kobe HPC Spring School12

机上で計算できる場合とおけば7-9, March 2022. Kobe HPC Spring School13

机上で計算できる場合を満足する2つのω（複素数）に対して,
とおきなる線形結合に対して、両項にを課せば、𝛼𝛼1 (−𝜛𝜛1𝑘𝑘−1 + 2𝜛𝜛1𝑘𝑘 − 𝜛𝜛1𝑘𝑘+1 ) + 𝛼𝛼2 ( − 𝜛𝜛2𝑘𝑘−1 + 2𝜛𝜛2𝑘𝑘 − 𝜛𝜛2𝑘𝑘+1 ) = 𝜆𝜆′ (𝛼𝛼1 𝜛𝜛1𝑘𝑘 + 𝛼𝛼2 𝜛𝜛2𝑘𝑘 )𝛼𝛼1 (−𝜛𝜛10 + 2𝜛𝜛11 − 𝜛𝜛12 ) + 𝛼𝛼2 ( − 𝜛𝜛20 + 2𝜛𝜛2 − 𝜛𝜛22 ) = 𝛼𝛼1 (2𝜛𝜛11 − 𝜛𝜛12 ) + 𝛼𝛼2 (2𝜛𝜛2 − 𝜛𝜛22 ) = 𝜆𝜆′ (𝛼𝛼1 𝜛𝜛1 + 𝛼𝛼2 𝜛𝜛2 )
𝛼𝛼1 (−𝜛𝜛1𝑁𝑁−1 + 2𝜛𝜛1𝑁𝑁 − 𝜛𝜛1𝑁𝑁+1 ) + 𝛼𝛼2 ( − 𝜛𝜛2𝑁𝑁−1 + 2𝜛𝜛2𝑁𝑁 − 𝜛𝜛2𝑁𝑁+1 )=𝛼𝛼1 (−𝜛𝜛1𝑁𝑁−1 + 𝜛𝜛1𝑁𝑁 ) + 𝛼𝛼2 ( − 𝜛𝜛2𝑁𝑁−1 + 𝜛𝜛2𝑁𝑁 ) = = 𝜆𝜆′ (𝛼𝛼1 𝜛𝜛1𝑁𝑁 + 𝛼𝛼2 𝜛𝜛2𝑁𝑁 )λが固有値, uが対応する固有ベクトルとなるより

7-9, March 2022. Kobe HPC Spring School

さらに

14

固有値とは
 工学における現象解析、システム解析に利用する

構造物の耐震振動解析

ネットワークなどの定常状態解析

 計算方法：線形代数の理論上は次の特性多項式を解けばよい



コンピュータ上で行列式の計算は難しい。多項式の求解は？ニュートン法でできるか？

7-9, March 2022. Kobe HPC Spring School15

ネットワーク解析 ある時刻にサイトをアクセスしている人が、次の時刻に別のサイトに移動する確率が、サイトとサイトの間で決まっているとしよう。
をサイトにいる人の数とする。

遷移確率に従って人が移動すれば定常状態(k∞)では

7-9, March 2022. Kobe HPC Spring School16

ネットワーク解析 実際に人が移動する状態を確認してみよう＊[100,0,0]の状態からスタート[100,0,0]->[40,0,60]->[34,18,48]->[31.6, 18, 50.4]->[31.36,18.72,49,92] -> [31.264,18.720,50.016]-> … -> [31.25, 18.75, 50.0]
定常的には100人のうちの半数の50人がサイト３を見ていることになる。
7-9, March 2022. Kobe HPC Spring School17固有「の」：特異値解析
 特徴量の主成分解析
ベクトルとして表現された複数データの主たる特徴



類似性の計算

平均的な顔の定義（例）固有「の」スキャンデータから切り出した１０個の「の」の最も平均的な特徴を示す「の」を計算7-9, March 2022. Kobe HPC Spring School18密行列向け固有値計算アルゴリズム

7-9, March 2022. Kobe HPC Spring School

19

べき乗法

7-9, March 2022. Kobe HPC Spring School20

（絶対値）最大固有値の計算方法
 べき乗法

初期ベクトルを選択



ベクトルに行列を乗じて正規化する



ベクトルが収束したら


ノルムが絶対値最大固有値対応する固有ベクトル7-9, March 2022. Kobe HPC Spring School

21

べき乗法
 行列

に
を乗じて、正規化し と置き直す。これを繰り返す。下に、実際にどのような計算が進行するかを示します。

3.15

3.163

3.164

7-9, March 2022. Kobe HPC Spring School

2.23

3.00

3.13

3.1642

3.1642

3.1642522

べき乗法 反復系統の解法の基本であるので、原理と解けるための条件をまと

める



条件：

の絶対固有値が全て異なるとする(最大だけで十分)

初期ベクトルは固有ベクトルが一次独立なので線形結合で書ける



ここで、次のように

7-9, March 2022. Kobe HPC Spring School

を更新していく

23

べき乗法
 ここで

は０でないと仮定すると、

．．．

とおく

の更新の方法から

7-9, March 2022. Kobe HPC Spring School

24

べき乗法
 ポイント

絶対値最大の固有値に対する固有ベクトルに収束

絶対値最大固有値に重複があると、一般に収束しない（複数のベクトルが張る空間内で変化し続けます）
絶対値最小固有値計算にも応用可能
が成り立つことから、逆行列の絶対値最大固有値は元の行列の絶対値最小固有値に対応

ある値(σ)に近い、固有値も計算可能

7-9, March 2022. Kobe HPC Spring School

25

ヤコビ法

7-9, March 2022. Kobe HPC Spring School

26

ヤコビ法
 行列

を

対称行列とする。

 “全て”の固有値・固有ベクトルの組を求めたい。

7-9, March 2022. Kobe HPC Spring School

27

ヤコビ法
 行列を逐次相似変換して対角行列に収束


非対角項の絶対値最大の２＊２小行列を取り出し



回転変換により対角行列に（アルゴリズムから）を定め、回転変換を左右より乗じる

7-9, March 2022. Kobe HPC Spring School

28

ヤコビ法
 行列を逐次相似変換して対角行列に収束


実際は、下の箇所で計算

回転行列は、

行列



この相似変換により、行列行列に拡大し、対角(青)に1をおくが変換されるのは、ピンク色線の箇所
は徐々に対角行列に近づく

数学的な収束の証明は後半に7-9, March 2022. Kobe HPC Spring School

29

ヤコビ法
 ポイント： 行列を逐次相似変換して対角行列に収束

対角要素以外の中から絶対値最大の要素を探し、その要素を基準とした２ｘ２行列を基本とした２ｘ２のGivens回転行列を求め作用させていく。

7-9, March 2022. Kobe HPC Spring School

30

ヤコビ法
 相似変換を繰り返すから

とおけば





の各列ベクトルは固有ベクトルでもある。

 つまり、ヤコビ法の手続きを繰り返した結果得る対格行列の要素が固有値であり、その位置に対応する
クトルである。

7-9, March 2022. Kobe HPC Spring Schoolの列ベクトルが対応する固有ベ31

ヤコビ法の原理
 はたして対角行列に収束しているのか？

を0に変換するとき、i,jを含む行と列に変更が加わる。

に対して、対角を除いた成分の２乗和を定めると→対角以外の部分の２乗和は単調減少 → 対角行列に収束する
7-9, March 2022. Kobe HPC Spring School32ハウスホルダー三重対角法＋α7-9, March 2022. Kobe HPC Spring School33ハウスホルダー三重対角化
 ハウスホルダー変換

ベクトル uを鏡線として、鏡面に反対の位置に移す変換：

対称行列とするようなハウスホルダー変換を定めることができる。
なるベクトルで、

7-9, March 2022. Kobe HPC Spring Schoolととればよい。

34

ハウスホルダー三重対角化
 ハウスホルダー変換

：

対称行列

とするようなハウスホルダー変換を定めることができる。
とブロック拡大表記した相似変換を定め、 に左右から作用に対しても同様の変換を再帰的に実施しすれば, 左下の成分は対角成分以下に非ゼロ要素をもつ行列、すなわち三重対角行列（対称）に最終的になる。
三重対角行列は要素数も少なく、固有値計算も容易。多くのアルゴリズムが存在する。
7-9, March 2022. Kobe HPC Spring School35

三重対角行列の固有値計算
 ハウスホルダー三重対角化後の行列Tと変形できるので、行列式を計算すると
つまり、Tの固有値はAの固有値と同じである。
 固有ベクトルは

Tの固有ベクトルに対して

7-9, March 2022. Kobe HPC Spring School

を作用（逆変換）すればよい

36

ＱＲ法
 ＱＲ分解

：

今,

上三角行列

,

,…

なる手続きを進めていくと、最終的にある行列に収束する性質を使い固有値・固有ベクトルを計算する。
ただし、は対称行列とする。
→

これより、

の固有値は

固有ベクトルはの各成分,の対応する列ベクトル

行列の更新手順ではのみ必要で、角の構造を保持するため計算が容易

7-9, March 2022. Kobe HPC Spring School

は常に三重対37

２分法
 三重対角行列の特性多項式の計算とおくと、一般形としてにより特性多項式を求められる。ゼロ点を所謂２分法により求めればよい。
ただし、nが大きいときはオーバフローの危険があり、ストゥルムの数え上げが有効。
（２分法）

なる区間(x,x’)に対して中点(x+x’)/2の符号により区間縮小を行う

7-9, March 2022. Kobe HPC Spring School

38

逆反復法
 何らかの方法でTの固有値近似が得られたとする。

、さらには固有ベクトルの近似

 べき乗法の箇所でも説明したように,に近い固有ベクトルはの絶対値最大固有値に対応する固有ベクトルである。

 従って,

により、固有ベクトルの精度を高める計算を行う。
ただし、本方法で求めた固有ベクトルは直交性の観点からはよくないので、直交化の手続きをしなくてはいけない

7-9, March 2022. Kobe HPC Spring School

39

分割統治法 三重対角行列を適当な摂動により以下のようにする

何らかの方法でとの固有値計算が為されたとする。それぞれの固有値と固有ベクトルを並べた行列(など) を用いてと問題を簡単化できる(相似変換なので固有値はもとのものと同じ、固有ベクトルは逆変換が必要)

7-9, March 2022. Kobe HPC Spring School40

分割統治法 簡単化された問題, ｛対角行列+1階摂動｝
適時変形を行っていけば次式を得る。

この有理方程式は、区間独立に解くことが容易である。

に解が存在することが分かっているので

 固有ベクトル

上記で求めたを用いて、以下の式より

7-9, March 2022. Kobe HPC Spring Schoolを計算する。最終的には正規化が必要

41

ＱＲ法もう一回
 先の説明でQR法を（対称）三重対角行列に適用したが、一般的な非対称行列に対しての有効な解法はQR法しかないのが現状である。

 非対称行列の場合はハウスホルダー変換を実施しても、三重対角化はできず、ヘッセンベルグ形式（下副対角成分以下が０の行列）に変形してから、前述のＱＲ法を行い上三角行列に収束する性質を利用して固有値を計算する。（複素固有対がある場合は2x2の対角ブロックが出現）
7-9, March 2022. Kobe HPC Spring School
42ハウスホルダー逆変換
 ハウスホルダー変換自身はそれ自身が逆変換でもあるので、作用させた逆順に作用していけばよい。

 近代的なアルゴリズムでは複数のハウスホルダー変換をまとめてブロック化する。

2つ以上の場合にも適用は可能である。なお、添字の順番などで形式が違うこともあるので注意。
この変形で、殆どの演算が行列と行列の積に置き換えられる。

7-9, March 2022. Kobe HPC Spring School43疎行列向けの反復解法
反復解法の多くは

Templates for the Solution of Algebraic
Eigenvalue Problems, A Practical Guide,
SIAMを参考にしています

7-9, March 2022. Kobe HPC Spring School44べき乗法から

7-9, March 2022. Kobe HPC Spring School

45

べき乗法
 疎行列では、記憶領域の問題やフォーマット変形操作が煩雑といった理由から、行列の変形操作を行わない。

 主に、行列とベクトルもしくはベクトル同士の積、近似もしくは写像し縮小した行列での操作が中心となる。

密行列でも扱ったが、べき乗法が計算の基本となる。

7-9, March 2022. Kobe HPC Spring School

46

べき乗法
 疎行列の場合、絶対値最大固有１つという計算はあまり行われず、最大値から数個、最小値から数個という場合が多い。

 複数ベクトルを用いて、べき乗法の原理から

7-9, March 2022. Kobe HPC Spring School47

ランチョス法の前に 先のべき乗法でが登場したが、これは直交基底
Xによって張られる部分空間の射影問題を扱っているといえる。

 今Ａの近似固有解に直交するように決めるに対して、その残差をその元となる部分空間

これは以下と同値です。
また,

が直交基底
となっており（

7-9, March 2022. Kobe HPC Spring Schoolで張られるのであれば、の線形結合）

48

Rayleigh/Ritz
 この様な、固有空間に対してガレルキン近似を入れて計算する方法を、

Rayleigh/Ritzの手法と呼びます。
分といえます。

に対する射影が本質部

 実際には、空間を張るm本の基底ベクトルから、射影した望するk(<m)個の近似固有値を取得する。

 次に、それに対する有ベクトルを得る。

の所の固有ベクトルにVを乗じて得られる近似固 このとき、上記の近似固有値をRitz値, 対応する近似固有ベクトルを

Ritzベクトルと呼ぶ。

7-9, March 2022. Kobe HPC Spring School49

ランチョス法 先の解法には, 固有ベクトルを近似する空間の張り方に自由度があった。ここで、クリロフ部分空間法によって生成される基底を使用してみる。

といった計算により, Vが逐次生成されていく。Ａの対称性を考慮すると

Ｔはαが対角、βが副対角に並ぶ三重対角行列である。rは反復の打ち切りで生じる項
7-9, March 2022. Kobe HPC Spring School50

ランチョス法（アルゴリズム）

7-9, March 2022. Kobe HPC Spring School51

ランチョス法 先の枠組みに当てはめれば, Tの固有値固有ベクトルがRitz値,

Ritzベクトルになっている。

 Ritz対によるAに対する残差を見てみると

 つまり, 反復を打ち切る際にでるのノルム係数
いかに、Ritzベクトルの精度はよっていることになる

7-9, March 2022. Kobe HPC Spring School

が十分に小さ52

ランチョス法（アルゴリズム）

7-9, March 2022. Kobe HPC Spring School53アーノルディ法 Aが対称行列（もしくはエルミート）であった条件を緩和して、非対称行列にする。基本的な枠組みは同様であり

 の形で固有ベクトルを形成する部分空間を作成していき、適当な精度が出たところで打ち切る。

 ただし、途中計算で生成される射影行列は、ヘッセンベルグ行列となる。

7-9, March 2022. Kobe HPC Spring School54アーノルディ法（アルゴリズム）

7-9, March 2022. Kobe HPC Spring School55

リスタート ランチョス法、アーノルディ法ともに 反復により部分空間を張る基底を作成し続けても、所望の固有ベクトルを十分近似できないことがある。

 その様な場合に、適切な部分空間基底を初期対として再度反復を開始するリスタートの技術が存在する。

 Thick Restart Lanzcosなどいくつかの手法があるが、今回は名前だけを紹介するのみとする。

（リスタートアルゴリズム）
Thick Restart, Implicit Restart, Explicit Restartなどなど

7-9, March 2022. Kobe HPC Spring School

56

Jacobi-Davidoson
 Krylov部分空間法とは異なる原理で、部分空間生成基底を生成する方法としてJD法が存在する。

Ritz対から生じる条件を課して基底を作る
実際には以下の、方程式を解き基底 t を作る

ここで、

は

の成分を排除するフィルタの役割をする。

7-9, March 2022. Kobe HPC Spring School57

JD法（アルゴリズム）

7-9, March 2022. Kobe HPC Spring School58

その他の解法
 対称行列の固有値問題はRayleigh商の極小問題

を局所最小化することでも計算可能  非線形関数の最小化問題
LOBPCGやCGR, MINRESなどの手法もあるが、今回は省略します。

7-9, March 2022. Kobe HPC Spring School59数値計算ライブラリ

7-9, March 2022. Kobe HPC Spring School60

数値計算ライブラリとは 高度なプログラミングを保証するための構成要素の一つです。
 非常に複雑な数学的機能、アルゴリズム、スキーム、データ処理などを完備なインターフェイスで提供します。
 例：連立一次方程式求解、高速フーリエ変換または離散フーリエ変換（FFT or DFT)、固有値計算、特異値計算、最小化、統計処
理など
 参照コードとして役立つことがあります 良い（悪い）プログラミングの例かもしれません***

Form C := alpha*A*B + beta*C.

50

DO 90 j = 1,nIF (beta.EQ.zero) THEN
DO 50 i = 1,mc(i,j) = zero
CONTINUEELSE IF (beta.NE.one) THEN
DO 60 i = 1,mc(i,j) = beta*c(i,j)

60

70
80
90

CONTINUE
END IF
DO 80 l = 1,ktemp = alpha*b(l,j)
DO 70 i = 1,mc(i,j) = c(i,j) + temp*a(i,l)
CONTINUE
CONTINUE
CONTINUE

http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html

7-9, March 2022. Kobe HPC Spring School61数値計算ライブラリとは 高度なプログラミングを保証するための構成要素の一つです。
 非常に複雑な数学的機能、アルゴリズム、スキーム、データ処理などを完備なインターフェイスで提供します。
 例：連立一次方程式求解、高速フーリエ変換または離散フーリエ変換（FFT or DFT)、固有値計算、特異値計算、最小化、統計処
理など
 参照コードとして役立つことがあります 良い（悪い）プログラミングの例かもしれません
 多くの場合、専門家により作成であり、あなたが作ったものよりも優れた計算速度と計算精度を提供してくれます 市販のライブラリ：より速く、より正確だが高価
 オープンソースライブラリ：高速で自由な利用が認められてるライセン
ス形態（商用ライブラリよりも高速な場合もある)
アプリケーションコードを実行する前にどれが皆さんに適しているかをチェックする必要があります。

7-9, March 2022. Kobe HPC Spring School

62

たとえば

 行列積がボトルネックになっているコードを高速化したい場合は、適切なBLAS (Basic Linear Algebra Subprograms) ライブラリの
DGEMMを使用してください (またはリンクしてください)。

Standard API for linear algebra kernels (case of (C)BLAS).

GEMM : Matrix-matrix multiplication()

AXPY: linear combination of 2 Vectors (
NRM2: Norm of a vector, etc. ()
Referencefor(i=0;i<N; i++) codes are available from netlib@UTK.
for(j=0;http://www.netlib.org/BLAS/
j<N; j++) {cblas_dgemm( CblasRowMajor,
)

t = 0.0;
CblasNoTrans, CblasNoTrans,for(k=0; k<N; k++)N, N, N, tCommercial:
Intel
MKL,
AMD
ACML
(free)+= a[i][k]*b[k][j];1.0, a, N, b, N, 0.0, c, N);c[i][j]=t;
 Open Source: ATLAS(@UTK), GotoBLAS(@TACC),} OpenBLAS for general purposed microprocessors
nVIDIA CUBLAS, AMD clMATH, MAGMABLAS(@UTK),KBLAS(@KAUST), ASPEN.K2(@RIKEN) : for GPGPU

25, September
7-9,
March 2022.
2020
Kobe
RIKEN
HPCCCS
Spring
HPCSchool
Summer School

63

たとえば

 行列積がボトルネックになっているコードを高速化したい場合は、適切なBLAS (Basic Linear Algebra Subprograms) ライブラリの
DGEMMを使用してください (またはリンクしてください)。

Standard API for linear algebra kernels (case of (C)BLAS).

GEMM : Matrix-matrix multiplication()

AXPY: linear combination of 2 Vectors (
NRM2: Norm of a vector, etc. ()
Referencefor(i=0;i<N; i++) codes are available from netlib@UTK.
for(j=0;http://www.netlib.org/BLAS/
j<N; j++) {cblas_dgemm( CblasRowMajor,
)

t = 0.0;
CblasNoTrans, CblasNoTrans,for(k=0; k<N; k++)N, N, N, tCommercial:
Intel
MKL,
AMD
ACML
(free)+= a[i][k]*b[k][j];1.0, a, N, b, N, 0.0, c, N);c[i][j]=t;
 Open Source: ATLAS(@UTK), GotoBLAS(@TACC),} OpenBLAS for general purposed microprocessors
nVIDIA CUBLAS, AMD clMATH, MAGMABLAS(@UTK),KBLAS(@KAUST), ASPEN.K2(@RIKEN) : for GPGPU

25, September
7-9,
March 2022.
2020
Kobe
RIKEN
HPCCCS
Spring
HPCSchool
Summer School

64

たとえば

 行列積がボトルネックになっているコードを高速化したい場合は、適切なBLAS (Basic Linear Algebra Subprograms) ライブラリの
DGEMMを使用してください (またはリンクしてください)。

Standard API for linear algebra kernels (case of (C)BLAS).
GEMM : Matrix-matrix multiplication()

AXPY: linear combination of 2 Vectors (
NRM2: Norm of a vector, etc. ()
Reference codes are available from netlib@UTK.
http://www.netlib.org/BLAS/
)

商用: Intel MKL, AMD ACML (free) オープンソース: ATLAS(@UTK), GotoBLAS(@TACC), OpenBLAS
for 汎用マルチコアプロセッサ向け
 nVIDIA CUBLAS, AMD clMATH, MAGMABLAS(@UTK),KBLAS(@KAUST), ASPEN.K2(@RIKEN) : GPGPU向け
25, September
7-9,
March 2022.
2020
Kobe
RIKEN
HPCCCS
Spring
HPCSchool
Summer School65

もっと紹介します
提案：より複雑な問題の場合は、以下を使用してください。


LAPACK (http://www.netlib.org/lapack/)

ScaLAPACK (http://www.netlib.org/scalapack/)Elemental (http://libelemental.org/)

EigenExa(http://www.aics.riken.jp/labs/lpnctrt/EigenExa_e.html)

Dense, General

Dense Eigenvalue



ELPA (http://elpa.rzg.mpg.de/)PETSc (http://www.mcs.anl.gov/petsc/)

Trillions (https://trilinos.org/)

ARPACK (http://www.caam.rice.edu/software/ARPACK/)FFTW (http://www.fftw.org/)

FFTE (http://www.ffte.jp/)



2decomp&FFT (http://www.2decomp.org/)MT, MTGP, dSFMT (http://www.math.sci.hiroshima-u.ac.jp/~mmat/MT/SFMT/index.html)
Random number

GMP big number librart(https://gmplib.org/)

QD pack, MPACK, and so on

7-9, March 2022. Kobe HPC Spring School

Sparse, General

Sparse, Eigenvalue
FFT

Multi-precision number
66

Rokkoでの実演習

7-9, March 2022. Kobe HPC Spring School

67

県立大Rokkoサーバでの実行時の注意事項
 環境変数LANGを適切に設定（抜けていることがあるので）% export LANG=C Intelコンパイラ等を利用しますので、ログインノード上で以下を実行して、モジュールを有効にしてください。

% module load intel intelmpi
petsc-intelmpi slepc-intelmpi gnuplot[guest31@rokko1 ~]$ module list
Currently Loaded Modulefiles:
1) intel/19.1.3
4) slepc-intelmpi/3.16.2
2) intelmpi/2020.4.304 5) gnuplot/5.2.7
3) petsc-intelmpi/3.16.2
7-9, March 2022.
Kobe HPC Spring School68ソースファイルとヒントは、~guest31/share/
SLEPc/  SLEPc sample code, a job-script, and Makefile Please ‘cd’ your work directory, then ‘cp -R’ thedirectory (‘%’ is a prompt, and do not type). You willsee
% cp -R ~guest31/share ./third_day% cd third_day
% ls
SLEPc

 他の例題ファイルも多数用意されてます(tests以下を変えて探してください)。

$PETSC_DIR/arch-linux-c-debug/tests/ksp/ksp/tutorials$SLEPC_DIR/arch-linux-c-debug/tests/eps/tutorials
7-9, March 2022. Kobe HPC Spring School

69

PETSc(SLEPC)
Official site on PETSc and SLEPC
Hands-on exercises

7-9, March 2022. Kobe HPC Spring School

70

PETSc/TAO
 Developed by Argonne National Lab. USA

Portable, Extensible Toolkit for Scientific Computation
 Toolkit for Advanced Optimizationhttps://www.mcs.anl.gov/petsc/
The current version of PETSc is 3.13; released March
29, 2020.
PETSc, pronounced PET-see (the S is silent), is a suiteof data structures and routines for the scalable (parallel)solution of scientific applications modeled by partial
differential equations. It supports MPI, and GPUs
through CUDA or OpenCL, as well as hybrid MPI-GPU
parallelism. PETSc (sometimes called PETSc/Tao) also
contains the Tao optimization software library.
(cf. PETSc/TAO homepage)
7-9, March 2022. Kobe HPC Spring School

71

SLEPc Significant plugin modules of PETSc for EVPdeveloped by Universitat Politecnica de Valencia, Spain

https://slepc.upv.es/
SLEPc is a software library for the solution of large
scale sparse eigenvalue problems on parallel computers.
It is an extension of PETSc and can be used for linear
eigenvalue problems in either standard or generalized
form, with real or complex arithmetic. It can also be
used for computing a partial SVD of a large, sparse,
rectangular matrix, and to solve nonlinear eigenvalue
problems (polynomial or general). Additionally, SLEPc
provides solvers for the computation of the action of amatrix function on a vector.
The current version of SLEPc is 3.13, released in March,
2020.
7-9, March 2022. Kobe HPC Spring School72

サンプルコードへのアクセス
https://slepc.upv.es/handson/handson1.html

リンク先URL が無効なので.
以下のディレクトリをコピーしてください
Other tutorial examples are on/work/gt65/t65021/share/SLEPc/work/gt65/t65021/share/PETSc

7-9, March 2022. Kobe HPC Spring School73コンパイルとリンク

makefile(一部のみ掲載)

CFLAGS =FFLAGS =CPPFLAGS =FPPFLAGS =LOCDIR = ${PWD}
EXAMPLESC = ex1.c
EXAMPLESF = ex1f.F
MANSEC = EPS
include ${SLEPC_DIR}/lib/slepc/conf/slepc_common

make コマンドを使います！！
make ex1でex1が生成されます！

7-9, March 2022. Kobe HPC Spring School74

プログラム実行（ジョブ投入）#!/bin/bash#PBS -q S#PBS -l select=1:ncpus=1#PBS -N slepc_exe#PBS -j oe
module load intel intelmpi petsc-intelmpi slepc-intelmpi
cd ${PBS_O_WORKDIR}
mpiexec.hydra -n ${NPROCS} ./ex1 -n 540 -m 180 -eps_type gd -eps_nev 4 eps_monitor -eps_view

 qsub run.shでジョブ投入します

7-9, March 2022. Kobe HPC Spring School75

C version
1-D Laplacian Eigenproblem, n=100
Number of iterations of the method: 19
Solution method: krylovschur
Number of requested eigenvalues: 1
Stopping condition: tol=1e-08, maxit=100
Number of converged eigenpairs: 2

C versionk||Ax-kx||/||kx||----------------- -----------------3.999033
4.02784e-09
3.996131
4.31174e-09
F90 version
1-D Laplacian Eigenproblem, n =100 (Fortran)
Number of iterations of the method: 19
Solution method: krylovschur
Number of requested eigenvalues: 1
Stopping condition: tol=1.0000E-08, maxit= 100
Number of converged eigenpairs: 2

F90 versionk||Ax-kx||/||kx||----------------- -----------------3.9990E+00
4.0278E-09
3.9961E+00
4.3117E-09

7-9, March 2022. Kobe HPC Spring School76

PETSc/SLEPCで遊んでみます チュートリアルのページを参考にして,$ ./ex1-n 400 -eps_nev 3 -eps_tol 1e-7
$ ./ex1 -n 400 -eps_nev 3 -eps_ncv 24$ ./ex1 -n 100 -eps_nev 4 -eps_type lanczos

1-D Laplacian Eigenproblem, n=400
Number of iterations of the method: 100
Solution method: krylovschur
Number of requested eigenvalues: 3
Stopping condition: tol=1e-07, maxit=100
Number of converged eigenpairs: 1
k||Ax-kx||/||kx||----------------- -----------------3.999939
9.48781e-08

1-D Laplacian Eigenproblem, n=400

1-D Laplacian Eigenproblem, n=100

Number of iterations of the method: 60
Solution method: krylovschur

Number of iterations of the method: 62
Solution method: lanczos

Number of requested eigenvalues: 3
Stopping condition: tol=1e-08, maxit=100
Number of converged eigenpairs: 5

Number of requested eigenvalues: 4
Stopping condition: tol=1e-08, maxit=100
Number of converged eigenpairs: 4k||Ax-kx||/||kx||----------------- -----------------3.999939
9.48494e-09
3.999754
7.19493e-09
3.999448
1.18552e-09
3.999018
6.43926e-10
3.998466
1.04213e-09

k||Ax-kx||/||kx||----------------- -----------------3.999033
9.95783e-09
3.996131
1.97435e-09
3.991299
9.15231e-09
3.984540
3.55339e-09

7-9, March 2022. Kobe HPC Spring School77サンプルコードを学習 (Pseudo code)
SlepcInitialize( PETSC_NULL_CHARACTER, ierr )MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, …., n, n, ierr )MatSetUp( A, ierr )(…. Calculation of matrix elements and others …..)
EPSCreate( PETSC_COMM_WORLD, eps, ierr )EPSSetOperators( eps, A, PETSC_NULL_OBJECT, ierr )EPSSetProblemType( eps, EPS_HEP, ierr )EPSSolve( eps, ierr )EPSGetEigenPair( eps, …… )EPSDestroy( eps, ierr )SlepcFinalize( ierr )7-9, March 2022. Kobe HPC Spring School

Modern style
1. Initialization
2. Create a Handle for data
3. Setup parameters
4. Create a Handle for solvers
5. Setup parameters
6. Kick the solver
7. Retrieve the results
8. Destroy handles
9. Finalize78行列データはどうやって設定するの? (C)PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

// Simple matrix format
Mat
A;
EPS
eps;
EPSType tname;
PetscReal tol, error, *values;

matrix ハンドラ作成

MatCreate( PETSC_COMM_WORLD, &A )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N )MatGetOwnershipRange(A, &Istart, &Iend )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES )MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY )MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY )7-9, March 2022. Kobe HPC Spring School79行列データはどうやって設定するの? (C)PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

// Simple matrix format
Mat
A;
EPS
eps;
EPSType tname;
PetscReal tol, error, *values;

行列サイズMxN

MatCreate( PETSC_COMM_WORLD, &A )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N )MatGetOwnershipRange(A, &Istart, &Iend )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES )MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY )MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY )7-9, March 2022. Kobe HPC Spring School80行列データはどうやって設定するの? (C)PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

// Simple matrix format
Mat
A;
EPS
eps;
EPSType tname;
PetscReal tol, error, *values;

mxn ブロック行列
valuesから行列Aの値がセットされます

MatCreate( PETSC_COMM_WORLD, &A )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N )MatGetOwnershipRange(A, &Istart, &Iend )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES )MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY )MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY )7-9, March 2022. Kobe HPC Spring School81行列データはどうやって設定するの? (C)PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

// Simple matrix format
Mat
A;
EPS
eps;
EPSType tname;
PetscReal tol, error, *values;

分散状態のデータをアセン
ブルして、行列データを適切な状態で分散します

MatCreate( PETSC_COMM_WORLD, &A )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N )MatGetOwnershipRange(A, &Istart, &Iend )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES )MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY )MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY )7-9, March 2022. Kobe HPC Spring School82行列データはどうやって設定するの?(in fortran90)
PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

! Simple matrix format
Mat
A
EPS
eps
EPSType tname
PetscReal tol, error, values(:)

行列サイズMxN

call MatCreate( PETSC_COMM_WORLD, A, ierr )call MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )call MatGetOwnershipRange(A, Istart, Iend, ierr )call MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)call MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)call MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)
7-9, March 2022. Kobe HPC Spring School83行列データはどうやって設定するの?(in fortran90)
PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

! Simple matrix format
Mat
A
EPS
eps
EPSType tname
PetscReal tol, error, values(:)

行列サイズMxN

call MatCreate( PETSC_COMM_WORLD, A, ierr )call MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )call MatGetOwnershipRange(A, Istart, Iend, ierr )call MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)call MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)call MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)
7-9, March 2022. Kobe HPC Spring School84行列データはどうやって設定するの?(in fortran90)
PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

! Simple matrix format
Mat
A
EPS
eps
EPSType tname
PetscReal tol, error, values(:)

mxn ブロック行列
valuesから行列Aの値がセットされます

call MatCreate( PETSC_COMM_WORLD, A, ierr )call MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )call MatGetOwnershipRange(A, Istart, Iend, ierr )call MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)call MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)call MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)
7-9, March 2022. Kobe HPC Spring School85行列データはどうやって設定するの?(in fortran90)
PETScは内部データフォーマットとインターフェースデータを柔軟に処理します。
PETScのデータ管理機構により、ユーザはメモリ上の実際の状態を通常見ることができません。行列Aはハンドラ変数で処理され、行列要素はクエリAPIなどを介してアクセスされます。

! Simple matrix format
Mat
A
EPS
eps
EPSType tname
PetscReal tol, error, values(:)

分散状態のデータをアセン
ブルして、行列データを適切な状態で分散します

call MatCreate( PETSC_COMM_WORLD, A, ierr )call MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )call MatGetOwnershipRange(A, Istart, Iend, ierr )call MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)call MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)call MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)
7-9, March 2022. Kobe HPC Spring School

86Hands-on time
必要なファイルをコピーしてください.
SLEPc/  SLEPc sample code, a job-script, and Makefile

7-9, March 2022. Kobe HPC Spring School87

少し複雑な問題にチャレンジ 2次元問題の固有値分析を：１次元の場合とおくと、２次元の場合は次のように
ただし は対応する大きさの単位行列

7-9, March 2022. Kobe HPC Spring School88少し複雑な問題にチャレンジ 行列のsetup方法が変わる
MatGetOwnershipRange(A,&Istart,&Iend);for (i=Istart;i<Iend;i++) {if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }MatSetValue(A,i,i,2.0,INSERT_VALUES);}
MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);【ポイント】Mat Aのサイズを(n*n)x(n*n)に変更する
また、
長さn*nのベクトルだが必要に応じて(n,n)の２次元配列と処理する。

MatGetOwnershipRange(A,&Istart,&Iend);for (i=Istart;i<Iend;i++) {PetscInt ib=i/n; PetscInt i0=i%n;
if ( ib>0 ) { MatSetValue(A,i,i-n,-1.0,INSERT_VALUES); }if ( ib<n-1 ) { MatSetValue(A,i,i+n,-1.0,INSERT_VALUES); }if ( i0>0 ) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }if ( i0<n-1 ) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }MatSetValue(A,i,i,4.0,INSERT_VALUES);}
MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
7-9, March 2022. Kobe HPC Spring School89Play with SLEPC
 変更したプログラム（出力部も変更済み）$ ./ex1-2d -n 10 -eps_nev 3 -eps_tol 1e-7
$ ./ex1-2d -n 10 -eps_nev 3 -eps_ncv 24

2-D Laplacian Eigenproblem, n=10
Number of iterations of the method: 4
Solution method: krylovschur
Number of requested eigenvalues: 3
Stopping condition: tol=1e-07, maxit=100
Number of converged eigenpairs: 4k||Ax-kx||/||kx||----------------- -----------------7.837972
5.4621e-09
7.601493
6.44101e-08
7.365014
1.2931e-08
2-D Laplacian Eigenproblem, n=107.228707
1.40955e-08$ ./ex1-2d -n 10 -eps_nev 4 -eps_type lanczos
2-D Laplacian Eigenproblem, n=10
Number of iterations of the method: 3
Solution method: krylovschur

Number of iterations of the method: 10
Solution method: lanczos

Number of requested eigenvalues: 3
Stopping condition: tol=1e-08, maxit=100
Number of converged eigenpairs: 5

Number of requested eigenvalues: 4
Stopping condition: tol=1e-08, maxit=100
Number of converged eigenpairs: 4k||Ax-kx||/||kx||----------------- -----------------7.837972
2.40527e-14
7.601493
1.39769e-13
7.365014
3.14186e-11
7.228707
6.57291e-11
6.992229
1.14673e-09

k||Ax-kx||/||kx||----------------- -----------------7.837972
1.60317e-09
7.601493
1.72603e-10
7.601493
3.30641e-09
7.365014
2.12916e-09

7-9, March 2022.

90結果を可視化してみよう 固有ベクトルをファイルに書き出してみる{}

下記プログラム行を、ex1.cの130行あたりfor(i=のループ内部に追加する(a.awk は書式をgnuplot方式にするための外部スクリプト：下部枠内）
char filename[256];
sprintf(filename,"eigenvector%03ld.txt",i);PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INDEX);
VecView(xr,viewer);
PetscViewerPopFormat(viewer);
PetscViewerDestroy(&viewer);if ( Istart==0 ) {char command[1024];
sprintf(command,"awk -v n=%ld -v out=v%03ld.data -f a.awk %s",n, i, filename);
BEGIN{ if(n==0)n=10; if(out=="")out="ev000.data"; }system(command);/^[0-9]*:/{}
gsub(/:/,""); N=$1;

7-9, March 2022. Kobe HPC Spring School}row = N % n; col = N / n;printf("%d %d %lf%c¥n", row, col, $2, (row==n-1 ? "¥n" : " ")) > out;91

結果を可視化してみよう
実行ジョブスクリプト: n=10として、0付近の固有値を5個計算# mpiexec.hydra –n ${NPROCS} ./ex1-2d -n 10 eps_nev 5 -eps_target 0（右図）Gnuplotでsplot# set pm3d
# splot ‘v004.data’.answer-2dというディレクトリに
２次元のプログラムの回答例があります。
プログラムとスクリプトを使って実行した可視化結果が本ページ右グラフです。
rokkoにログインする際にssh –XによるXポート
フォワーディングしておく必要があります。
7-9, March 2022. Kobe HPC Spring School92

7-9, March 2022. Kobe HPC Spring School

93

まとめ
 代表的な並列数値計算ライブラリ


ScaLAPACK

EigenExa

PETSc

他にも多数存在。逐次版でLAPACK, BLASなど
 PETSc and SLEPc on Rokko

疎行列だけでなくで手軽に並列計算可能

アルゴリズムも選択できる



手続きが決まっているので、サンプルコードを参考にして必要な部分を追加拡張していく。
実行時の引数やプログラムコード内で設定変更可能

7-9, March 2022. Kobe HPC Spring School94