熱伝導方程式\[\frac{\partial u}{\partial t}=\Delta u\]を, 初期値\[u(t=0,x,y)=f(x,y)\]境界条件\[u(y=0)=1, u(y=1)=u(x=0)=u(x=1)=0\]で時間発展を追うことを考えましょう.

陽的解法

差分式は, 時間\(t\)の区間幅を\(\delta\)とすれば, \[\frac{u^{(n+1)}_{i,j}-u^{(n)}_{i,j}}{\delta}=\frac{u^{(n)}_{i-1,j}+u^{(n)}_{i+1,j}+u^{(n)}_{i,j-1}+u^{(n)}_{i,j+1}-4u^{(n)}_{i,j}}{\Delta^2}\]すなわち\[u^{(n+1)}_{i,j}=(1-4c)u^{(n)}_{i,j}+c\left(u^{(n)}_{i-1,j}+u^{(n)}_{i+1,j}+u^{(n)}_{i,j-1}+u^{(n)}_{i,j+1}\right),\quad c=\frac{\delta}{\Delta^2}\]ということになりますね.

• プログラムを組んで解を求めてみましょう. 解答例はこちら.

ですが差分式をよくみると,\[1-4c<0\]の場合，すなわち\[\frac{\delta}{\Delta^2}\geq\frac{1}{4}\]で，係数が負のものが現れるので, 雰囲気やばそうです. 実際, Fourier解析すると, 安定解が保証できないことを確かめることができるでしょう．

• 数学はまあ講義で勉強することにして，ここでは，cが大きすぎるとうまくいかないことを計算して確かめましょう．

これも, Courant-Friedrichs-Levyの条件と呼ばれます．

名前はともあれ，これは制限が厳しいです．空間解像度が高い計算がしたいな，と思って, まあ最初は\(\Delta=0.01\)くらいか？と思ったら, 時間ステップが\(\delta<0.000025\), つまり\(t=5\)まで求めようとしたら, 20万ステップ必要になります. これでは計算に時間がかかりすぎるですね！これでは4K, 8Kレベルの解像度\(\Delta=0.0001\)とかは夢のまた夢，ってことになります．

陰的解法

陽的解法では，計算が遅くて死にそうだと言うことがわかりました．そこで，陰的解法(インプリシット・スキーム)を考えましょう．ようするにこうします：\[ u^{(n+1)}_{i,j} = u^{(n)}_{i,j} +\frac{\delta}{\Delta^2}\left(u^{(n+1)}_{i-1,j}+u^{(n+1)}_{i+1,j}+u^{(n+1)}_{i,j-1}+u^{(n+1)}_{i,j+1}-4u^{(n+1)}_{i,j}\right)\]つまり\[u^{(n+1)}_{i,j}-\alpha\left(u^{(n+1)}_{i-1,j}+u^{(n+1)}_{i+1,j}+u^{(n+1)}_{i,j-1}+u^{(n+1)}_{i,j+1}\right) =\beta  u^{(n)}_{i,j}, c=\frac{\delta}{\Delta^2},\ \alpha=\frac{c}{1+4c},\ \beta=\frac{1}{1+4c}\]を解けば良い. これを行列で書いてみましょう.

\[ \begin{bmatrix}1 & -\alpha & & -\alpha\\ -\alpha & 1 & -\alpha & & \ddots\\ & -\alpha & 1 & \ddots & & \ddots\\ -\alpha & & -\alpha & \ddots & \ddots & & \ddots\\ & \ddots & & \ddots & \ddots & \ddots & & \ddots\\ & & \ddots & & \ddots & \ddots & \ddots & & \ddots\\ & & & \ddots & & \ddots & \ddots & \ddots & & \ddots\\ & & & & {\color{purple}-\alpha} & & {\color{magenta}-\alpha} & {\color{red}1} & {\color{orange}-\alpha} & & {\color{teal}-\alpha}\\ & & & & & \ddots & & -\alpha & \ddots & -\alpha & & \ddots\\ & & & & & & \ddots & & & \ddots & -\alpha & & -\alpha\\ & & & & & & & \ddots & & \ddots & \ddots & \ddots\\ & & & & & & & & \ddots & & \ddots & \ddots & -\alpha\\ & & & & & & & & & -\alpha & & -\alpha & 1 \end{bmatrix}\begin{bmatrix}u_{1,1}\\ \vdots\\ \vdots\\ \vdots\\ {\color{purple}u_{i-1,j}}\\ \vdots\\ {\color{magenta}u}_{{\color{magenta}i,j-1}}\\ {\color{brown}u_{i,j}}\\ {\color{orange}u_{i,j+1}}\\ \vdots\\ {\color{teal}u_{i+1,j}}\\ \vdots\\ u_{N-1,N-1} \end{bmatrix} \]\[=\begin{bmatrix}u_{1,1}^{(n)}\\ \vdots\\ \vdots\\ \vdots\\ {\color{purple}\beta u_{i-1,j}^{(n)}}\\ \vdots\\ {\color{magenta}\beta u^{(n)}}_{{\color{magenta}i,j-1}}\\ {\color{brown}\beta u_{i,j}^{(n)}}\\ {\color{orange}\beta u_{i,j+1}}\\ \vdots\\ {\color{teal}\beta u_{i+1,j}^{(n)}}\\ \vdots\\ u_{N-1,N-1}^{(n)} \end{bmatrix}+\begin{bmatrix}+\alpha u_{0,1}+\alpha_{1,0}\\ +\alpha u_{0,2}\\ \vdots\\ \vdots\\ 0\\ \vdots\\ \vdots\\ \vdots\\ 0\\ +\alpha u_{0,N-1}+\alpha u_{1,,N}\\ \vdots\\ \vdots\\ +\alpha u_{N-2,N-1}+\alpha u_{N-1,N-2} \end{bmatrix} \]

非斉次項は, 前時間ステップの値と, 境界条件の値から構成されているので注意．

あとは，これを解けば良い．とはいえ

• \(\Delta=0.01\)でも, 行列は結構でかい. \(10^4\times10^4\)ですからね．\(\Delta=0.0001\)であると,  行列サイズは\(10^8\times10^8\), ってことはメモリーサイズは\[\frac{8\times10^8\times10^8}{10^9}=8\times10^7\ \text{GB}\]ですね！富嶽スーパーコンピュータでも，これではできない．
• しかし, よくみると行列には5本しか筋がないので, 必要なメモリーは，本当は\[\frac{8\times10^8\times5}{10^9}=4\ \text{GB}\]しかないんですよ．これ, あんたのノートPCに入るかもしれないくらいじゃないか．

どうやら，疎行列をどう取り扱うのか，学ばねばならないようです．

SPARSE Matrix

疎行列というのは, まばらにしか非ゼロ要素が存在しない行列です．

まあでも, 数値シミュレーションに出現する疎行列では，対角要素はゼロではない，という条件がつくのがほとんどですので, ここでもそれを仮定しましょう．

疎行列の困ったちゃんなところは, 「疎行列を解こうとして未知数を消去(つまりLU分解)していったら，密行列になっちゃうかもよ」という点です．

\[\begin{bmatrix}\ 1\ & -1 & -1 & -1 & -1 & -1\\ 1 & 1\\ 1 & & \text{1}\\ 1 & & & 1\\ 1 & & & & 1\\ 1 & & & & & 1 \end{bmatrix}\to\begin{bmatrix}\ 1\ & -1 & -1 & -1 & -1 & -1\\ 0 & 2 & 1 & 1 & 1 & 1\\ 0 & 1 & 2 & 1 & 1 & 1\\ 0 & 1 & 1 & 2 & 1 & 1\\ 0 & 1 & 1 & 1 & 2 & 1\\ 0 & 1 & 1 & 1 & 1 & 2 \end{bmatrix}\]

疎行列が溢れた・・・このままじゃ俺の計算は全滅じゃあ〜〜 
落ち着いて．あれは・・・Eigenに向かって飛べ！　　　　　

疎行列自身の要素数（＝計算機パワー）は十分小さくても，連立方程式を解こうとしてLU分解を始めると，要素数が爆発して計算機パワーを凌駕してしまうわけね．

しかし，まあ数値シミュレーションに出現する疎行列では，ブロック状に0行列が出現することが多い:

この場合，理屈で言って「ここゼロだし」ということが分かることも多いわけで，そこを回避するようにプログラムを組むわけです．

上の例だと, A1-B1ブロックとC4-A4ブロックの計算では，オレンジ部分はずっと0のままである

だが，それは，問題によってプログラムが特注品になる，ということでもあります．だから，頑張ってね．ただし,

• Eigenには, 疎行列を取り扱うためのライブラリーがあるみたい

なので，手はあるのだろうと思います．まあでも，どれくらい爆発するものか経験しておいた方が良いかもなので，ここでは，自家製教育用疎行列ライブラリー libSPARSE を使ってみましょう.

ここで作成したmyArrayクラスとlibSPARSEを使った解答例は，こちら.

\(\Delta=0.01\)でもノートPCでなんとか実行できます．\(\delta=0.1\)でも発散しないけど, ビデオをきれいに撮るために\(\delta=0.001\)で実行したのが下のやつです．

XCodeで実行中に使用メモリーサイズを調べると，LU分解している時に使用メモリーが微妙に爆発していることがわかります．どれくらい爆発するかは，問題によって異なるでしょう．このケースでlibSPARSEであると，数倍に増えるみたいですね．

インプリシットスキームは\(\delta\to\infty\)でも発散しないので，極限をとってみましょう. すると\[4u^{(n+1)}_{i,j}-\left(u^{(n+1)}_{i-1,j}+u^{(n+1)}_{i+1,j}+u^{(n+1)}_{i,j-1}+u^{(n+1)}_{i,j+1}\right) =0\]になって, 定常のLaplace方程式を解いていることになります．

GITレポジトリーを使う

さて，いよいよ手元のノートPCでは計算量が限界に近づいてまいりました．今後は

• 研究室の計算機クラスター
• 諸大学のスーパー計算機クラスター

を利用することになっていくでしょう．そうなると，プログラムやデータの通信受け渡しが問題になります．とりあえず問題になるのはプログラムのほうですね．これを簡単に行うのがGITレポジトリです．

GITとは何か，ここのページを読んで, とりあえずBitbucketのアカウントの取得を行っておいてください．そしてここの【パスワードなしGIT】の【公開鍵登録】で，あなた自身のノートPCの公開鍵（ここで作ったやつ）を登録しておいてください．

ビルドシステムCMake

XCodeでは，プログラムの再生には再生ボタンをクリック，となっていました．計算機クラスターでは, GUI(グラフィカル・ユーザーインターフェース)がありませんので，再生ボタンがありません．その対策として，ビルドツールを使うことになります．

ここでは, CMakeを使ってみましょう．まずMacにおいて，XCodeの再生ボタンではなく，【ターミナル】で再生が行えるようにします．

まずは再生ボタンソフトCMakeをここに従ってインストール．ほとんどの人は不要と思うが, pkgconfigもインストールして損はない．

まずXCodeのプロジェクトのフォルダーをFinderで開いてください：

 

こんな感じですかね．ここに, 以下のCMakeLists.txt というファイルをエディターで作成してください．

てゆうか，なんでかそれが面倒くさいんですよね．

• こんな対策をしている人がいます．
• 私は下の方で【新規ターミナル】を開いて touch CMakeLists.txt とコマンドを入力してしまいます．
  • フォルダーに【新規ターミナル】がでてこない場合はここを参照

で，こうなる

Finderにでてきたら，クリックでもして編集します．作成するファイルは，例えばこんなものです:

cmake_minimum_required(VERSION 3.9)
set(CMAKE_BUILD_TYPE Release CACHE STRING "None Debug Relase ...")
#ここには，あなたのプロジェクト名を指定
project(many_programs)
#言語標準ば設定
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_C_STANDARD 99)
set(CMAKE_C_STANDARD_REQUIRED ON)
#パッケージば検索
find_package(Boost 1.70.0 REQUIRED COMPONENTS
  date_time coroutine date_time filesystem
  graph log_setup program_options random
  regex serialization system thread timer
)
include_directories(${Boost_INCLUDE_DIRS})
link_directories(${Boost_LIBRARY_DIRS})
#あんたのターゲット名(複数可能)を記述していく
#   多分ターゲット名フォルダにmain.cppプログラムがあるので, それを追加. cppをリストアップする
add_executable(SPARSE
  SPARSE/main.cpp
)
target_link_libraries(SPARSE
  ${Boost_LIBRARIES}
)
#他にもあるならどんどん追記

これで準備完了です．

Macビルド

ターミナルでプロジェクトの場所にいるとしましょう. そこで次のように入力します：

$ pwd
/Users/sugimoto/Documents/..../many_programs  ← 念のため，確認しましょうね．
$ mkdir Darwin            ← folder名は他でも良いけど
$ cd Darwin
$ cmake ..                ← cmake 空白 点点 ですよ！
-- The C compiler identification is AppleClang 12.0.5.12050022
-- The CXX compiler identification is AppleClang 12.0.5.12050022
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /Library/Developer/CommandLineTools/usr/bin/cc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /Library/Developer/CommandLineTools/usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Found Boost: /usr/local/lib/cmake/Boost-1.71.0/BoostConfig.cmake (found suitable version "1.71.0", min


imum required is "1.70.0") found components: date_time coroutine date_time filesystem graph log_setup pro


gram_options random regex serialization system thread timer
-- Found PkgConfig: /usr/local/bin/pkg-config (found version "0.29.1")
-- Found GSL: /usr/local/include (found version "2.6")
-- Found MPI_C: /usr/local/lib/libmpi.dylib (found version "3.1")
-- Found MPI_CXX: /usr/local/lib/libmpi.dylib (found version "3.1")
-- Found MPI: TRUE (found version "3.1")
-- Configuring done
-- Generating done
-- Build files have been written to: /Users/sugimoto/Documents/Source/XCode/many_programs/Darwin

最後に, Build files have been written to とか，うまくいってそうな雰囲気が漂っていることを確認してください．標準的なライブラリがインストールされていると, 勝手に検索してしまうので，メッセージには個人差があります．

では再生しましょう．

$ make SPARSE    ← あなたのターゲット名を指定
[ 50%] Building CXX object CMakeFiles/SPARSE.dir/SPARSE/main.cpp.o
[100%] Linking CXX executable SPARSE
[100%] Built target SPARSE

これで「ビルド」までが完了です．フォルダを ls コマンドで見ると, ターゲット名のファイルができています：

$ pwd
/Users/sugimoto/Documents/..../many_programs/Darwin ← 念のため，確認しましょうね．
$ ls -l
total 312
-rw-r--r--   1 sugimoto  staff   27523  6  4 11:45 CMakeCache.txt
drwxr-xr-x  14 sugimoto  staff     448  6  4 11:45 CMakeFiles
-rw-r--r--   1 sugimoto  staff    5489  6  4 11:45 Makefile
-rwxr-xr-x   1 sugimoto  staff  118579  6  4 11:45 SPARSE      ← こいつ！これが実行ファイル
-rw-r--r--   1 sugimoto  staff    1584  6  4 11:45 cmake_install.cmake

では，再生（プログラムの実行）は，あなたが行います：

$ ./SPARSE   ← 実行ファイルの前に ./   (てんわる)をつけるのがミソ
nt= 1 err= 0.558864
nt= 2 err= 0.172615
nt= 3 err= 0.0586968
nt= 4 err= 0.0197901
nt= 5 err= 0.00665965
nt= 6 err= 0.00224012
nt= 7 err= 0.000753444
nt= 8 err= 0.000253408
nt= 9 err= 8.52287e-05
nt= 10 err= 2.8665e-05

データファイルは，プログラムであなたが指定した場所に作成されますので，Darwinフォルダにできるとは限らないので注意.

このように, Macでも，【ターミナル】があれば, プログラムを再生できます．というわけで，計算機クラスターでも，完全に同じ手順を踏みます．それには CMakeLists.txt が必要です．そこで, XCodeにそれを登録しておきましょう．

【プロジェクト】で【Add Files】

 【CMakeLists.txt】選択

で, CMakeLists.txt が出現するので，【Source Control】【Add CMakeLists.txt】

GITレポジトリとXCode合体

ここに従って, あなたのXCodeのプロジェクトとGITレポジトリを合体させましょう．合体が完了した場合，以下の利点が得られます．

• あなたのノートPCが消え失せても，あなたの仕事はWEBに残っている．
  • 新しい機体をゲットしたら数分で作業を続行することが可能である．
• お定まりのイベント「修士論文提出期限3日前にPCが動かなくなる」でも，問題はない
  • 修論間際で作業量アップすると，大体壊れるんだよね・・・
  • 大切なプログラムはWEBにあるので, PCはゴミ箱にポイで問題ない
• そもそも，ノートPCで作業をする必要がないので，計算機クラスターで作業をすることもできる
• 今日は気分でベッドの上でノートPCで入力したなら，それもよい

というわけで，計算機クラスターだろうがスーパーコンピュータだろうが，もはや大した違いは無くなったわけです．

計算機クラスターを使う

計算機クラスターにログインし，XCodeで書いてきたプログラムをダウンロードしましょう．BitBucketで【クローンの作成】でコマンドをゲット：

ほんだば，行ってみましょう. 計算機クラスターにログインからね．

arm64:MyMac sugimoto$ ssh sun0ext   ← VPNの人は ユーザー名@sun1 とかだな
fucker135sugimoto@forward.kuins.kyoto-u.ac.jp's password:うにゃにゃにゃ
Last login: Thu Jun  3 13:45:29 2021 from 10.249.229.55
[sugimoto@sun0 ~]$ mkdir example   ← フォルダー名は各自で選定してね
[sugimoto@sun0 ~]$ cd example/
[sugimoto@sun0 example]$ git clone git@bitbucket.org:sugimoto605/many_programs.git
Cloning into 'many_programs'...
remote: Counting objects: 45, done.
remote: Compressing objects: 100% (41/41), done.
remote: Total 45 (delta 5), reused 0 (delta 0)
Receiving objects: 100% (45/45), 26.43 KiB | 1.32 MiB/s, done.
Resolving deltas: 100% (5/5), done.
[sugimoto@sun0 example]$

ダウンロードできたので, ビルドして再生ね:

[sugimoto@sun0 example]$ cd many_programs/    ← あなたのプロジェクト名
[sugimoto@sun0 many_programs]$ mkdir Linux_gcc
[sugimoto@sun0 many_programs]$ cd Linux_gcc/
[sugimoto@sun0 Linux_gcc]$ cmake ..
-- The C compiler identification is GNU 8.4.1
-- The CXX compiler identification is GNU 8.4.1
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /usr/bin/cc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Found Boost: /usr/local/lib/cmake/Boost-1.71.0/BoostConfig.cmake (found suitable version "1.71.0", m


inimum required is "1.70.0") found components: date_time coroutine date_time filesystem graph log_setup


 program_options random regex serialization system thread timer
-- Found PkgConfig: /usr/bin/pkg-config (found version "1.4.2")
-- Found GSL: /usr/include (found version "2.5")
-- Found MPI_C: /usr/local/lib/libmpi.so (found version "3.1")
-- Found MPI_CXX: /usr/local/lib/libmpi.so (found version "3.1")
-- Found MPI: TRUE (found version "3.1")
-- Configuring done
-- Generating done
-- Build files have been written to: /net/sun0/sugimoto/example/many_programs/Linux_gcc

最後に, Build files have been written to とか，うまくいってそうな雰囲気が漂っていることを確認してください．

[sugimoto@sun0 Linux_gcc]$ make SPARSE  ← あんたのターゲット名だよ
Scanning dependencies of target SPARSE
[ 50%] Building CXX object CMakeFiles/SPARSE.dir/SPARSE/main.cpp.o
[100%] Linking CXX executable SPARSE
[100%] Built target SPARSE
[sugimoto@sun0 Linux_gcc]$ ls -l　　← 一応確認
合計 132
-rw-rw-r-- 1 sugimoto fluid 27569  6月  4 12:38 CMakeCache.txt
drwxrwxr-x 6 sugimoto fluid   320  6月  4 12:40 CMakeFiles
-rw-rw-r-- 1 sugimoto fluid  5345  6月  4 12:38 Makefile
-rwxrwxr-x 1 sugimoto fluid 93184  6月  4 12:40 SPARSE    ← こいつ！これが実行ファイル
-rw-rw-r-- 1 sugimoto fluid  1667  6月  4 12:38 cmake_install.cmake
[sugimoto@sun0 Linux_gcc]$ ./SPARSE
nt= 1 err= 0.558864
nt= 2 err= 0.172615
nt= 3 err= 0.0586968
nt= 4 err= 0.0197901
nt= 5 err= 0.00665965
nt= 6 err= 0.00224012
nt= 7 err= 0.000753444
nt= 8 err= 0.000253408
nt= 9 err= 8.52287e-05
nt= 10 err= 2.8665e-05

これでLinuxでも動くことができました．今のところ，ファイルサーバーの貧相なCPU(8個くらい)で計算を行っています．長大な（数十時間）かかるような計算は，ジョブ管理システムの利用を検討し，計算専用ムキムキマシン(数十個以上のCPU)の利用を検討してくださいね．ファイルサーバーは磁石ディスクですが，計算サーバーの新しげなのは全部SSDなので，ファイル読み書きも桁違いです．スーパーコンピュータはさらに二桁ほど違います．

もちろん，ただ単に計算機CPUをたくさん占有しても，計算が速くなるわけではありません．

• 並列計算機コードを書く
• 並列計算機コンパイルをする

必要があります．もはや初心者向けではありませんね．こちらのOpenMPを参照してください．