sakeisnomitaiの日記

プログラミングで調べたことを自分なりにまとめるだけの人

Open MPを用いて並列計算 バグの出現編

前回の記事で並列計算を試みたことを書いた。簡単な並列計算はこれで十分できそうなのだが、私が研究で使っているプログラムがなぜか上手く動かない。

sakeisnomitai.hatenablog.com

 

今回はこのプログラムの並列化をするのが目標となる。実際のプログラムは研究内容と関わるため割愛。7重ループを用いたFORTRAN77のプログラムである。

 

目標:7重ループの最も外側のループのみを並列化。多分このループは独立しているはず。

 

1. エラーその1

前回のようにとりあえず、プログラムの初めに

program ○○

       use omp_lib

       implicit double precision(a-h, o-z)

のようにする。(implicit noneを使えというツッコミはとりあえずなしで。)

また、並列化させたい部分に

!$omp parallel

!$omp do

      do 50 j = 1, nj

というように元のプログラムに追記していく。

 

これで実行するとSegmentation faultという始末。。。

 

解決策1 : スタックサイズの変更

メモリのスタックサイズが足りないためにsegmentation faultとなることは並列計算をするときにはよくあるらしい。とりあえずの解決法として、ターミナルで

ulimit -s unlimited

とすることで解決できる。unlimitedの部分を具体的な数値で指定することもできる。

現在の設定を見るには、

ulimit -a

とすれば良い。

 

→結果: 計算は回ったものの、数値が発散してしまった。

 

解決策2 : プログラムの中身の検討

結局はここに戻ってくる。まずは並列化させたいループの初めでデバッグ用のwrite文を追加

 

do 50  jy = 1,njy

 ys=(jy-0.5)*isy

※ここでisy=1

 (途中数行省略)

write(6, *) 'debug' , omp_get_thread_num(), jy, ys, isy

 

こうすると出力は...

 

 

(3列目) - 0.5 = (4列目)になるはずなのに、そうはなっていないことがわかった。

これは、各スレッドごとに変数が共有されていて、上書きが生じていることが原因だと推測できる。

 

・つまり、各スレッドごとに独立した変数を持つように設定すれば良い

 

調べたところこのようなものが見つかった。NAG万歳。

www.nag-j.co.jp

 

このページによると、並列計算の開始時に変数の属性を宣言すれば良い。

つまり

 

!$ omp parallel private(jy,ys)

!$ omp do

 

というように書けば問題のys, jy を各スレッドで独立に扱うことができる。(他のスレッドと共有させたい場合には private の代わりに shared を用いる)

 

変数の属性を指定してあげた場合、以下のような出力を得られる。

 

jy, ys についてはうまく設定できたことがわかる。(但し計算結果は発散)

 

他の変数についてもprivate属性を設定すれば上手く計算できそうだ。

 

 

まとめ

並列計算を導入するときにそのままomp do文を付け足すだけだと、変数が各スレッドで共有されてしまって予期せぬ計算結果が出てしまうことがあるので、それを独立させることが必要!

 

 

 

 

 

 

 

Macで並列計算できるようにしたい(FORTRAN) (記事更新 2/13)

Macでちょっと時間のかかる数値計算をすると待ち時間がめんどい。。。3時間かかる計算は嫌だ。。。

 

ということで並列化に手を出してみます。

どうやらOpen MPかOpen MPIを使えばいいみたいです。何が違うの?

 

・Open MP

   どうやら1台のPC内の複数スレッドで計算できるようにするものらしい

・Open MPI 

 複数のプロセッサ間で通信して並列計算できるみたい

www.cc.u-tokyo.ac.jp

 

ここのサイトがわかりやすそう

 

 

とりあえずOpen MPを使えるようにしよう!

まずは、私の環境について。

M1 Macbook Pro 2021

OSはVentura

今回はFortranを想定(一応私はFortran90使いです)

 

Open MPをインストールする

 

特にインストールするものなんて無かった!!(gccがあればOK)

brew install libomp 

が必要でした

 

 

コンパイルオプションで 

  -fopenmp 

をつけてあげれば大丈夫

 

Open MPプログラムの作り方

まず冒頭部分に

use omp_lib

どOpen MPを使うことを宣言する。

 

並列領域の初めに

!$omp parallel

並列領域の終わりに

!$omp end parallel

と書くことで、並列化する場所だけを明示する。

 

例えば、あるループの部分を並列計算させたい場合には、

 

!$omp parallel

!$omp do

do i = 1, n

    (doループの処理)

end do

!$omp end do

!$omp end parallel

 

と書けば良い。"omp do"を書かないと全てのスレッドで全てのループを実行してしまうらしい。また、1スレッドのみで処理を行いたい部分については"!$ omp single"を用いれば良い。

 

 

プログラムの実行

export OMP_NUM_THREADS=4

などとしてスレッド数を設定(bashの場合, この場合4スレッドで計算)

 

これをやればできる。。。はず!

 

 

 

参考

http://www.eccse.kobe-u.ac.jp/assets/images/simulation_school/kobe-hpc-summer-basic-2018/KHPC2018-openmp.pdf

https://www.cc.kyushu-u.ac.jp/scp/doc/users/lecture/2018/openmp-nov2018.pdf

OpenMP入門: 初めてのOpenMPプログラム - Hello World in OpenMP

 

 

追記(20230213)

なぜか1行目のprogram ○○の部分でsegmentation faultしてしまったが、bash

ulimit -s unlimited

として、スタックサイズの制限をなくすことでとりあえず動くようになった。が、数値が発散しているので更なるお勉強が必要です。。。次の記事はこれかな。

 

 

 

 

 

MacOS Venturaでgfortranが動かなくなった話

Mac初心者なもので、うっかり何も考えずにMac OSをVenturaにアップデートしてしまった。するとなぜかgfortranが使えない...

 

具体的には、以下のようなエラーが出た。

 

ld: -rpath can only be used when targeting Mac OS X 10.5 or later

collect2: error: ld returned 1 exit status

 

とりあえずgfortran --versionでバージョン確認すると、

gfortran: warning: could not understand version '13.00.00'

と帰ってくる。???

 

ここでようやくMac OSの更新が悪さをしていることに気づいたので、諸々のアップデートをする。(Xcodeは最新版になっていたので省略)

 

1. homebrewの更新

brew upgrade  でOK

 

そしたらちゃんとコンパイルできるようになった。

誰かに一助となれば幸いです。

 

[追記]

pythonも併用していますがpython, pipはアップグレードする必要がありました。

また、numpyなどのモジュールもpipでインストールし直さないとダメでした。

OSのアップデートは計画的に(自戒)。

Macでlapackを使う

連立一次方程式を解くのにlapackという数値計算用のフリーのライブラリがあることを知ったので、ちょっと使ってみた。環境はM1Macです。

 

1.インストール

 →homebrewでできる

brew install lapack  でOK

 

インストールすると...

For compilers to find lapack you may need to set:

  export LDFLAGS="-L/opt/homebrew/opt/lapack/lib"

  export CPPFLAGS="-I/opt/homebrew/opt/lapack/include"

と言われる。この通りexportしてあげれば良いみたい。

 

2.連立一次方程式を解いてみる

倍精度で解くには DGESVというサブルーチンを使う

n次方程式:Ax = b

 

call DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)

で呼び出せる。

N, LDA, LDBはnで大丈夫ぽい

NRHSはBの列数 1でOK

A:n×nの係数行列 B: 1×nのベクトル(既知)

IPIV: integer,dimension(n)

Info: integer

 

このサブルーチンを呼び出すとbに解が格納される。

write(6,*) b 

としてしまえば解を出力できる。

 

コンパイル

gfortran ~.f90 -llapack

とすれば出来た。

 

詳しく知りたければ

www.nag-j.co.jp

とか

 

sxauroratsubasa.sakura.ne.jp

とか

qiita.com

とかを見てみると良さそうです。(参考にさせていただきました)

はじめに

勉強日記始めました。勉強したことを吐き出すためだけの日記です。

 

Fortran90での数値計算について、またはPythonでのお絵描きについて困ったことをまとめて投稿していくつもりです

 

世界のどこかにいる同じことでお困りの人に届けー!!!!