数値計算ライブラリの使用方法「固有値計算法を中心にして」
Kobe HPC Spring School ２０１９
今村俊幸 理化学研究所 計算科学研究センター
Toshiyuki Imamura, RIKEN Center for Computational Science
2019/3/13～ 15

本日の講義（１） 固有値（線形代数）と応用問題

振動問題
 ネットワーク定常問題
 固有値計算アルゴリズム
 密行列






べき乗法
ヤコビ法
ハウスホルダー三重対角＋分割統治法＋逆変換

疎行列




ランチョス法
ヤコビ・デビッドソン法
その他 固有値計算ソフトウェア紹介

ScaLAPACK EigenExa
 PETScはじめに表記法について
 ベクトル：

小文字で表記（’→’記号はつけない）。
m行(縦方向) 行列：大文字で表記する。
m行(縦方向), n列(横方向)の様に, ベクトルを並べた記法使用も 数体：

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

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

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

 行列のノルム：

 行列式：
 条件数：

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

縦１本線表記表記法について 固有値固有ベクトル

固有値をギリシャ記号のラムダ( ) もしくはシータ( )固有ベクトルはで注意。

もしくは を用いて表現するが、適時記号を変えるの



固有値と対応する固有ベクトルを「固有対」と呼ぶ



固有値の記号に対応した大文字の行列は、固有値を対角に並べてできる対角行列とする



また、固有ベクトルを並べてできるベクトルも同様に同じ文字の大文字で表記することがある。

固有値(線形代数)と応用問題

主に線形代数の復習から

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

構造物の耐震振動解析

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

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



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

固有値計算の例
 柱の変形

境界条件を考慮して解は

多自由度系の振動
 支配方程式系



変数分離法で解ける



行列・ベクトルで整理するとの形をしている

机上で計算できる場合とおけば

机上で計算できる場合を満足する2つのω（複素数）に対して,とおきなる線形結合に対して、両項にを課せば、λが固有値, uが対応する固有ベクトルとなるよりさらに

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

構造物の耐震振動解析

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

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



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

ネットワーク解析 ある時刻にサイトをアクセスしている人が、次の時刻に別のサイトに移動する確率が、サイトとサイトの間で決まっているとしよう。
をサイトにいる人の数とする。
遷移確率に従って人が移動すれば定常状態(k∞)ではネットワーク解析 実際に人が移動する状態を確認してみよう＊[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人がサイト３を見ていることになる。

固有「の」：特異値解析
 特徴量の主成分解析
ベクトルとして表現された複数データの主たる特徴



類似性の計算

平均的な顔の定義（例）固有「の」スキャンデータから切り出した１０個の「の」の最も平均的な特徴を示す「の」を計算

密行列向け固有値計算アルゴリズム

べき乗法

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

初期ベクトルを選択



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



ベクトルが収束したら


ノルムが絶対値最大固有値対応する固有ベクトルべき乗法
 行列に を乗じて、正規化し と置き直す。これを繰り返す。下に、実際にどのような計算が進行するかを示します。

3.15

3.163

3.164

2.23

3.00

3.13

3.1642

3.1642

3.16425

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

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

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



ここで、次のようにを更新していくべき乗法 ここでは０でないと仮定すると、 の更新の方法から

．．．

とおく

べき乗法
 ポイント

絶対値最大の固有値



絶対値最大固有値に重複があると、一般に収束しない（複数のベクトルが張る空間内で変化し続けます）



絶対値最小固有値計算にも応用可能


に対する固有ベクトルに収束が成り立つことから、逆行列の絶対値最大固有値は元の行列の絶対値最小固有値に対応

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

ヤコビ法

ヤコビ法
 行列を対称行列とする。

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

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

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



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

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

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

回転行列は、

行列



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

数学的な収束の証明は後半に

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

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

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

とおけばの各列ベクトルは固有ベクトルでもある。

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

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

に対して、対角を除いた成分の２乗和を定めると→対角以外の部分の２乗和は単調減少 → 対角行列に収束する

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

ベクトル uを鏡線として、鏡面に反対の位置に移す変換：対称行列とするようなハウスホルダー変換を定めることができる。
なるベクトルで、ととればよい。

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

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

Tの固有ベクトルに対してを作用（逆変換）すればよい

ＱＲ法
 ＱＲ分解：上三角行列,今,,…

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

これより、の固有値は固有ベクトルはの各成分,の対応する列ベクトル

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

は常に三重対角の構

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

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

逆反復法
 何らかの方法でTの固有値近似、さらには固有ベクトルの近似 が得られたとする。
 べき乗法の箇所でも説明したように,に近い固有ベクトルはの絶対値最大固有値に対応する固有ベクトルである。

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

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

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

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

この有理方程式は、区間独立に解くことが容易である。に解が存在することが分かっているので

 固有ベクトル

上記で求めたを用いて、以下の式よりを計算する。最終的には正規化が必要ＱＲ法もう一回
 先の説明でQR法を（対称）三重対角行列に適用したが、一般的な非対称行列に対しての有効な解法はQR法しかないのが現状である。
 非対称行列の場合はハウスホルダー変換を実施しても、三重対角化はできず、ヘッセンベルグ形式（下副対角成分以下が０の行列）に変形してから、前述のＱＲ法を行い上三角行列に収束する性質を利用して固有値を計算する。（複素固有対がある場合は2x2の対角ブロックが出現）







ハウスホルダー逆変換
 ハウスホルダー変換自身はそれ自身が逆変換でもあるので、作用させた逆順に作用していけばよい。

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

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

疎行列向けの反復解法
反復解法の多くは
Templates for the Solution of Algebraic
Eigenvalue Problems, A Practical Guide,SIAMを参考にしていますべき乗法から

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

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

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

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

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

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

 今Ａの近似固有解に対して、その残差をその元となる部分空間に直交するように決める
これは以下と同値です。
また,が直交基底となっており（で張られるのであれば、の線形結合）

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

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

に対する射影が本質部分といえます。

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

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

ランチョス法 先の解法には, 固有ベクトルを近似する空間の張り方に自由度があった。

ここで、クリロフ部分空間法によって生成される基底を使用してみる。

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

Ｔはαが対角、βが副対角に並ぶ三重対角行列である。rは反復の打ち切りで生じる項

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

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

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

 つまり, 反復を打ち切る際にでるのノルム係数に、Ritzベクトルの精度はよっていることになるが十分に小さいかランチョス法（アルゴリズム）

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

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

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

アーノルディ法（アルゴリズム）

リスタート ランチョス法、アーノルディ法ともに 反復により部分空間を張る基底を作成し続けても、所望の固有ベクトルを十分近似できないことがある。
 その様な場合に、適切な部分空間基底を初期対として再度反復を開始するリスタートの技術が存在する。

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

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

Jacobi-Davidoson
 Krylov部分空間法とは異なる原理で、部分空間生成基底を生成する方法

としてJD法が存在する。

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

ここで、はの成分を排除するフィルタの役割をする。

JD法（アルゴリズム）

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

を局所最小化することでも計算可能  非線形関数の最小化問題
LOBPCGやCGR, MINRESなどの手法もあるが、今回は省略します。
時間があれば、LOBPCGの追加資料を用意しておきます。

数値計算ライブラリ密行列ソルバー
 有力なものを列挙します。
1.

EISPACK (1970年代の古いものです)2.

LAPACK (逐次計算の業界標準)

3.

PLASMA (multiCore環境向け）

4.

MAGMA(GPU/MIC向けのもの）

5.

ScaLAPACK (分並列計算での標準)

6.

Elemental

7.

ELPA

8.

EigenExa

疎行列ソルバー
 有力なものを列挙します
1.

ARPACK

2.

PARPACK

3.

PETSc (SLEPC)

4.

Trilinos5.

Z-Paresライブラリ使用法（概要編）
実習編資料は別途配布

ScaLAPACKの利用方法

ScaLAPACKサイトの例題より
 http://www.netlib.org/scalapack/examples/

ｓample_pdsyev_call.fをダウンロードする.

How to use ScaLAPACK
 Link方法, Intel compiler+MPIの環境で% mpiifort –o exe ｓample_pdsyev_call.f –lscalapack –lmkl_rt
 実行:#!/bin/bash#PJM -L "rscgrp=school"#PJM -L "node=2x2"#PJM -L "elapse=00:05:00"#PJM –mpi “proc=4”#PJM -jexport OMP_NUM_THREADS=16
mpiexec ./exe

Let’s learn the sample code
 行列生成関数: PDLAMODHILB
 固有値計算関数: PDSYEV 結果出力: PDLAPRNT

他に、初期化(BLACS_XXX)や終了(BLACS_EXIT)、行列データを扱うためのdescriptorの宣言(DESCINIT)などが必要。
BLACS: 通信回りの関数系、通常利用者からはプロセスの２次元配置の仕方などの管理系と思えばよい
PDSYEV: ＱＲ法による固有値計算ルーチン
他に、PDSYEVX, PDSYEVD, PDSYEVRなど存在する行列データはプロセス間で「２次元ブロック分割」されており、行列要素ごとに格納されるプロセスが決まっている。

行列設定関数を読む
 PDLAMODHILBを読んで分かるように

並列用の特別な変更は, 代入操作を関数呼び出しPDELSETにしている点。
PDELSETはもし、呼び出しプロセスがオーナーであれば指定された値をローカルメモリ上の配列データにストアする

SUBROUTINE PDLAMODHILB( N, A, IA, JA, DESCA, INFO )DO 20 J = 1, N
DO 10 I = 1, N
IF( I.EQ.J ) THENCALL PDELSET( A, I, J, DESCA, ¥( DBLE( N-I+1 ) ) / DBLE( N )+ONE / ( DBLE( I+J )-ONE ) )ELSECALL PDELSET( A, I, J, DESCA, ONE / ( DBLE( I+J )-ONE ) )
ENDIF
END

やってみよう！
 PDLAMODHILBを改良して、自由な対称行列を入力として固有値計算をしてみよう。

 また、固定化されている行列次元(N), プロセス数（NPROW, NPCOL）を変えて、計算時間の変化を見てみよう！

 更に、余裕のある人は固有値求解関数をpdsyevdなどに変更してどうなるかを調べてみよう。

EigenExaの使用方法

EigenExaのダウンロード
http://www.aics.riken.jp/labs/lpnctrt/EigenExa.html

最新版をダウンロードしてくださHow to use EigenExa
 Build & Link方法 Intel compiler+MPIの環境で自身でリンクする% mpiifort –o exe foo.f –lEigenExa –lscalapack –lmkl_rt -lm 実行:#!/bin/bash#PJM -L "rscgrp=school"#PJM -L "node=2x2"#PJM -L "elapse=00:05:00"#PJM –mpi “proc=4”#PJM -jexport OMP_NUM_THREADS=16
mpiexec ./eigenexa_benchmark

Let’s learn the sample code
 行列生成関数: mat_set
 固有値計算関数: eigen_sx

ScaLAPACK同様に、初期化(eigen_init)や終了(eigen_free)必要。
EigenExaは2次元サイクリックサイクリック分割（具体的にはCOL,ROW共にNB=1固定の場合）の範囲ではScaLAPACKとコンパチであるので、相互の関数を利用し合うことができる。
従って、NB=1で行列のdescriptorが宣言されているので、PDLMODHILBで作成した配列をEigenExaに渡して解くこともできる。

やってみよう！
 行列次元(N), プロセス数（NPROW, NPCOL）を変えて、計算時間の変化を見てみよう！

 更に、余裕のある人はScaLAPACKのサンプルコード内のPDLAMODHILBを組み込んで、様々な行列の固有値求解をEigenExaにさせてみよう。

PETSc(SLEPC)の使い方
SLEPCの公式Hands-on exercisesサイトの題材を利用する

第一部は簡単な導入のみです
第二部ではPETScの使い方をもう少し深めます。

http://slepc.upv.es/handson/handson1.html

ex1.cをダウンロードしてください

ex1ｆ.Ｆをダウンロードしてくだ（スクロールした下の方にありますコンパイル&リンク
Makefile

ARCH = arch-fujitsu-sparc64fx-opt
PETSC_DIR = /opt/aics-ss/petsc-3.3-p5
SLEPC_DIR = /opt/aics-ss/slepc-3.3-p4
INCPATH= -I$(PETSC_DIR)/include -I$(PETSC_DIR)/$(ARCH)/include ¥-I$(SLEPC_DIR)/include -I$(SLEPC_DIR)/$(ARCH)/include
LDFLAGS= -L$(SLEPC_DIR)/$(ARCH)/lib -lslepc ¥-L$(PETSC_DIR)/$(ARCH)/lib -lpetsc ¥-SSL2BLAMP
all: ex1 ex1f
ex1: ex1.ompifccpx -o ex1 ex1.o $(LDFLAGS)
ex1f: ex1f.o
mpifrtpx -o ex1f ex1f.o $(LDFLAGS)
ex1.o: ex1.c
mpifccpx -c ex1.c -Xg $(INCPATH)
ex1f.o: ex1f.F
mpifrtpx -c ex1f.F $(INCPATH)
clean:¥rm ex1 ex1.o ex1f ex1f.omakeコマンドでコンパイルする ex1, ex1fが作成される

プログラムの実行 実行#!/bin/bash#PJM -L "rscgrp=school"#PJM -L "node=2x2"#PJM -L "elapse=00:05:00"#PJM –mpi “proc=4”#PJM -jexport OMP_NUM_THREADS=16
echo ‘C-version’mpiexec ./ex1 –n 100
echo ‘F90-version’
mpiexec ./ex1f –n 100

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バージョンの結果

k||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
k||Ax-kx||/||kx||----------------- -----------------3.9990E+00
4.0278E-09
3.9961E+00
4.3117E-09

F90バージョンの結果

Play with SLEPC
1-D Laplacian Eigenproblem, n=400

 チュートリアルページから$ ./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

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: 4

k||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

Learn the sample code
SlepcInitialize( PETSC_NULL_CHARACTER, ierr )MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, …., n, n, ierr )MatSetUp( A, ierr )(行列やベクトルデータの宣言とデータ設定)
ESPCreate( PETSC_COMM_WORLD, eps, ierr )ESPSetOperators( eps, A, PETSC_NULL_OBJECT, ierr )EPSSetProblemType( eps, EPS_HEP, ierr )EPSSolve( eps, ierr )EPSGetEigenPair( eps, …… )EPSDestroy( eps, ierr )SlepcFinalize( ierr )How to setup a matrix?
 PETScは内部データを柔軟に制御している。PETScが管理するため、使用者からは実態は見えない。行列ハンドラ変数 A を使ってアクセスする。内部フォーマットはいくつか存在する。
! Simple matrix format
MatAEPSeps
EPSType tname
PetscReal tol, error, values(:)
MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )MatGetOwnershipRange(A, Istart, Iend, ierr )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)
MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)
MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)

How to setup a matrix?
 PETScは内部データを柔軟に制御している。PETScが管理するため、使用者からは実態は見えない。行列ハンドラ変数 A を使ってアクセスする。内部フォーマットはいくつか存在する。
! Simple matrix format
MatAEPSeps
EPSType tname
PetscReal tol, error, values(:)

行列データの形成

MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )MatGetOwnershipRange(A, Istart, Iend, ierr )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)
MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)
MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)

How to setup a matrix?
 PETScは内部データを柔軟に制御している。PETScが管理するため、使用者からは実態は見えない。行列ハンドラ変数 A を使ってアクセスする。内部フォーマットはいくつか存在する。
! Simple matrix format
MatAEPSeps
EPSType tname
PetscReal tol, error, values(:)

行列サイズの指定
グローバルサイズ
MxN

MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )MatGetOwnershipRange(A, Istart, Iend, ierr )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)
MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)
MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)

How to setup a matrix?
 PETScは内部データを柔軟に制御している。PETScが管理するため、使用者からは実態は見えない。行列ハンドラ変数 A を使ってアクセスする。内部フォーマットはいくつか存在する。
! Simple matrix format
MatAEPSeps
EPSType tname
PetscReal tol, error, values(:)

mxnのブロック行列に対して配列valuesをセットする

MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )MatGetOwnershipRange(A, Istart, Iend, ierr )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)
MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)
MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)

How to setup a matrix?
 PETScは内部データを柔軟に制御している。PETScが管理するため、使用者からは実態は見えない。行列ハンドラ変数 A を使ってアクセスする。内部フォーマットはいくつか存在する。
! Simple matrix format
Mat
AEPS
eps
EPSType tname
PetscReal tol, error, values(:)

セットされた配列
データをアセンブルする

MatCreate( PETSC_COMM_WORLD, A, ierr )MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, M, N, ierr )MatGetOwnershipRange(A, Istart, Iend, ierr )MatSetValues( A, m, idxm, n, idxn, values, INSERT_VALUES|ADD_VALUES, ierr)
MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY, ierr)
MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY, ierr)２部でやってみよう！
 少し難しくなりますが、長方領域板振動を表現する方程式、重調和関数の方程式(E: ヤング率, h: 板厚, ρ: 密度, μ: ポアソン比)
変数分離法により次の固有方程式を得る

やってみよう！
 固定境界条件(ディリクレ条件)のもと差分化してやってみよう！
 テンソル積の表示を使って

 簡単化のため正方形状として、Nx=Ny=m, hx=hy=1とする.

やってみよう！
 この様に作られる行列Aの固有値と固有ベクトルを計算して、得られた固有値の大きい方から順に選んで対応する固有ベクトルをm×mに並べ換えたデータとして可視化するとどうなるだろうか？
 可視化にはgnuplotのplot3dを使ってみましょう。