Fortranのコンパイルオプション（fortran77/90で有効）

デバッグオプション

上記のオプションをつけて実行すると詳しいエラーメッセージが表示されて計算が止まる。

最適化オプション

計算のちょっとしたコツ

ライブラリの使用

ライブラリとしてMath Kernel Library (MKL)をインストールしている。ここでは線形代数計算(Lapack)に対してMKLを使用する例を紹介する。バージョンが少し古いがMKLの日本語マニュアルのリンクをここに貼っておく。マニュアルに従ってプログラム内でLapackのサブルーチンをコールし、コンパイルする際にMKLへリンクを指定することでLapackのルーチンを使うことができる。具体的には 
ifort -L /opt/intel/compilers_and_libraries/linux/mkl/lib/intel64/ -mkl hoge.f90
のようにコンパイルオプションを指定するだけでLapackのルーチンを使うことができる。MKLには乱数生成やフーリエ変換もサポートしている（詳しい使い方はマニュアルを参照）。

並列化マニュアル

並列計算を用いる方法として、自動並列(OpenMP)と手動並列(OpenMPI)がある。それぞれの長所と短所は以下のとおり。

長所 短所
自動並列(OpenMP) プログラムちょっと書き換えてコンパイルオプションを付け加えれば簡単に並列化できる。並列効率が良くなるよう自動で調整してくれる PC間をまたがった並列計算（ノード間並列）ができず、1台のPCの複数のコアを用いた並列計算（ノード内並列）のみ可能。全て自動なので、細かい並列処理が難しい。
手動並列(OpenMPI) ノード内並列・ノード間並列の両方が可能である。並列の仕方について細かく指定できるので、OpenMPでは並列化できなかった部分を並列化できたりもする。メモリを分散させるための並列計算が可能である。 並列化の仕組みとMPI言語を知る必要があり、プログラムを組むのに大きな労力がかかる。並列効率を意識してプログラムを書く必要があり、下手なプログラムを書くと全く速くならなかったりする。

OpenMP並列

ここではよく使う2つの例を挙げておく。 OpenMPを用いたプログラムをコンパイルするにはコンパイルオプション-openmpを用いる。
ifort -openmp hoge.f90
また、実際に並列計算をさせるにはOMP_NUM_THREADSという環境変数にコア数を与えなければいけない。ジョブスクリプトの例としては 
#!/bin/bash
#PBS -q batch
#PBS -e error.log
#PBS -o out.log
#PBS -l nodes=1:ppn=4
#PBS -l walltime=00:10:00:00
cd $PBS_O_WORKDIR
export OMP_NUM_THREADS=4
./out.exe
となり、4行目のppnと7行目の環境変数の部分に並列コア数を指定する。また$PBS_NUM_PPNを用いて 
#!/bin/bash
#PBS -q batch
#PBS -e error.log
#PBS -o out.log
#PBS -l nodes=1:ppn=4
#PBS -l walltime=00:10:00:00
cd $PBS_O_WORKDIR
export OMP_NUM_THREADS=$PBS_NUM_PPN
./out.exe
としてもよい。

OpenMPI並列

ここではよく使う部分だけ記述する。詳しくは理研のマニュアルにて。

MPI並列ではすべてのコアが、同一のプログラムを実行することになる。最も簡単なOpenMPIプログラムの例を示す。 
program main
  implicit none
  include '/opt/openmpi/include/mpif.h'
  integer :: ierr,nprocs,myrank
!
  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,nprocs,ierr)
  call mpi_comm_rank(mpi_comm_world,myrank,ierr)
!
  write(*,*) myrank,nprocs
!
  call mpi_finalize(ierr)
  stop
end program main
10行目でmyrankとnprocsを標準出力に書き出している。このプログラムを4並列で実行すると 
0 4
1 4
2 4
3 4
となる（表示される行の順番はバラバラかもしれない）。これは4つのコアが独立に上記のプログラムを実行し、それぞれのコアが10行目に来た時点でmyrankとnprocsを出力した結果である。MPI並列計算をさせるために書かれたプログラムをコンパイルするにはifortの代わりにmpif90（FORTRAN77で書いているのであればmpif77を用いる）を用いる。MPI並列計算をさせるためのジョブスクリプトの例として 
#!/bin/bash
#PBS -q batch
#PBS -e error.log
#PBS -o out.log
#PBS -l nodes=4:ppn=1
#PBS -l walltime=00:10:00:00
cd $PBS_O_WORKDIR
export OMP_NUM_THREADS=1
mpirun ./out.exe
となる。実行にはmpirunを用いること。ここでnodesとppnの違いに注意。ppnはノード内でどれだけ並列させるかというのを指定し、nodesはノード間でどれだけ並列させるかというのを指定する。例えばnodes=1:ppn=4とすると、mauiは現時点で最も速いノードを探し出し、そのノードに対して並列数4のノード内並列計算をさせる。逆にnodes=4:ppn=1とすると1並列ずつノードにジョブを割り振りながらノードの速度を随時測定する。コアの割り振りが決まった時点で、すべてのコアが同一ノードに割り振られていればノード内並列になる（つまりnodes=1:ppn=4としたときとほぼ同じになる）し、複数のノードに割り振られていれば、ノード間並列となる。全並列数=nodes×ppnとなり、例えばnodes=2:ppn=4とするとノード内4並列を各2ノードで行うような8並列の計算となる（nodes=2に対してmauiが同じノードに計算を割り当てれば、ノード内8並列になる）。ノード間並列ではLANケーブルを介したデータ通信を行うため、データ通信が多ければ多いほど並列効率は悪くなる。

以下に具体的なプログラム例を示す。

完全に独立なdo文の並列化

いわゆるバカパラと呼ばれる並列計算で、つまり並列化しなくても、ジョブを個別に複数投入すれば代用できる類の並列化である。例えば異なるパラメーターで全く同じ計算をさせたい場合に 
program main
  implicit none
  include '/opt/openmpi/include/mpif.h'
  integer :: ierr,nprocs,myrank
!
  integer :: it
  real(8) :: t
  他の変数を定義
!
  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,nprocs,ierr)
  call mpi_comm_rank(mpi_comm_world,myrank,ierr)
!
  do it=myrank,100,nprocs
    t=dble(it)/100.d0
    計算
  enddo
!
  call mpi_finalize(ierr)
  stop
end program main
これはパラメータ0≤t≤1においてtを0.01刻みに動かして計算させる際に、パラメータの動かし方を並列化させたという最も簡単な並列化である。この並列計算において注意したい点は、101通りのtの計算が各コアごとに均等に割り振られるため、特にノード間並列をさせた場合に速いコアは先に終わってずっと待っており、遅いコアがずっと計算することになる。これを回避するために次のような改良が考えられる。 
program main
  implicit none
  include '/opt/openmpi/include/mpif.h'
  integer :: ierr,nprocs,myrank
  character(len=20) :: crank
!
  integer :: it
  real(8) :: t
  他の変数を定義
!
  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,nprocs,ierr)
  call mpi_comm_rank(mpi_comm_world,myrank,ierr)
!
  it=myrank
  if(myrank==0) then
    open(10,file='next.dat',form='formatted')
    write(10,'(i5)') nprocs
    close(10)
  endif
  write(crank,*) myrank
  call system('sleep '//crank)
!
  do
    t=dble(it)/100.d0
    計算
    open(10,file='next.dat',form='formatted',status='old')
    read(10,'(i5)') it
    close(10)
    if(it>100) exit
    open(10,file='next.dat',form='formatted')
    write(10,'(i5)') it+1
    close(10)
  enddo
!
  call mpi_finalize(ierr)
  stop
end program main
これは各tにおいてdo文の計算が終わったコアから次のtの計算を順番に行っていくためのプログラムである。ここでnext.datというファイルを作成しているが、ここには次に計算すべきitの値を書き出しておく。do文の一番最後まで来たコアは、next.datを読み込んで次のitの値を知り、next.datにit+1を上書きするという寸法である。非線形のシステムはnfsのファイル共有を随時更新するような設定になっているので問題ないが、nfsのファイル共有にタイムラグをもたせているようなシステムの場合、この方法は使えない（next.datの更新が遅れてすでに終わっているitの計算を再度別のコアが行う可能性があるため）。next.datへの書き込みおよび読み込みに対して複数のコアが集中するのを防ぐためにdo文の直前で(myrank)秒だけ計算を止め、コアごとに計算のタイミングを少しずらしている。

Monte Carlo計算などでアンサンブルを稼ぎたい場合に、各アンサンブルの計算を並列化させるという並列計算が考えられる。この場合は平均値等を計算するために、最後にそれぞれのコアで計算された値を足し合わせる必要がある。前の例を一歩進めて書き直せば 
program main
  implicit none
  include '/opt/openmpi/include/mpif.h'
  integer :: ierr,nprocs,myrank
  character(len=20) :: crank
!
  integer :: it
  real(8) :: t,s,st
  他の変数を定義
!
  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,nprocs,ierr)
  call mpi_comm_rank(mpi_comm_world,myrank,ierr)
!
  it=myrank
  if(myrank==0) then
    open(10,file='next.dat',form='formatted')
    write(10,'(i5)') nprocs
    close(10)
  endif
  write(crank,*) myrank
  call system('sleep '//crank)
!
  s=0.d0
  do
    計算
    s=s+st
!
    open(10,file='next.dat',form='formatted',status='old')
    read(10,'(i5)') it
    close(10)
    if(it>100) exit
    open(10,file='next.dat',form='formatted')
    write(10,'(i5)') it+1
    close(10)
  enddo
  mpi_allreduce(mpi_in_place,s,1,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
!
  if(myrank==0) then
    open(11,file='average.dat',form='formatted')
    write(11,'(e26.16e3)') s/101.d0
    close(11)
  endif
!
  call mpi_finalize(ierr)
  stop
end program main
これは101個のアンサンブルを並列化して計算させ、アンサンブル毎の値stを各コアで足し合わせる(=s)。最後にコア毎のsをすべて足し合わせる必要があるが、ここでデータ通信を行う。そのためのコマンドがmpi_allreduceであり、これは 
mpi_allreduce(通信したい変数,結果を格納する変数,通信したいデータ数,データの型,通信する際に行う演算,mpiコミュニケータ,エラーメッセージ)
の形で用いる。通信したい変数と結果を格納する変数を同じにしたい場合、通信したい変数の部分をmpi_in_placeとする。こうすると各ノードのsの和を計算して各ノードのsに格納する。また通信するデータは配列であっても良い。その場合は通信したい配列の一番最初のアドレスを指定し、通信したいデータ数を指定する。データの型としてmpi_integer,mpi_double_precision,mpi_double_complexなどがある。また上の例ではデータ通信の際にそれぞれのコアのデータを足し合わせるように指定したが、その場合にはmpi_sumとなる。その他にもmpi_productのように掛けたりする演算など色々指定できる。また、配列を指定した場合は配列の要素ごとに和をとったり積をとったりする。詳しくはマニュアルを参照。 基本的に通信を行うためのコマンドを指定すると、すべてのコアがそのコマンドの部分に来た時点で（つまり速いコアはそこで待っている）実行されるが、意図的にコアの同期（コアの足並みをそろえる）を行いたい場合には
mpi_barrier(mpi_comm_world,ierr)
とする。

独立ではないdo文の並列化：差分法

差分法であれば、普段ならOpenMPで全く問題ないが、乱数生成等との兼ね合いや、メモリの節約等でどうしてもOpenMPIを使いたい場合もあるので載せておく。また、差分法はMPI並列コードの格好の練習にもなる。今、2次元格子Nx X Nyを考え、Nyを並列化することを考える。また並列数をMとしたときにNy/Mは割り切れるとしておく。 
program main
  implicit none
  include '/opt/openmpi/include/mpif.h'
  integer :: ierr,nprocs,myrank
  integer :: istatus(mpi_status_size)
  integer :: kup,kdown
!
  integer,parameter :: nx=100,ny=100
  integer,parameter :: nxm=nx-1,nym=ny-1
  integer :: ix,iy
  integer :: nyi,nyf,nyr
  real(8),allocatable :: f(:,:),fd(:,:)
!
  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,nprocs,ierr)
  call mpi_comm_rank(mpi_comm_world,myrank,ierr)
!
  kup=myrank+1
  kdown=myrank-1
  if(myrank==0) kdown=mpi_proc_null
  if(myrank==nprocs-1) kup=mpi_proc_null
!
  nyr=ny/nprocs
  nyi=myrank*nyr
  nyf=(myrank+1)*nyr-1
!
  allocate(f(-1:nx,nyi-1:nyf+1),fd(0:nxm,nyi:nyf))
!
  do iy=nyi,nyf
    do ix=0,nxm
      f(ix,iy)=fに値を代入
    enddo
  enddo
!
  do iy=nyi,nyf
    f(-1,iy)=0.d0
    f(nx,iy)=0.d0
  enddo
  if(myrank==0) then
    do ix=0,nxm
      f(ix,-1)=0.d0
    enddo
  else if(myrank==nprocs-1) then
    do ix=0,nxm
      f(ix,ny)=0.d0
    enddo
  endif
!
  call mpi_sendrecv(f(0,nyf),nx,mpi_double_precision,kup,1,f(0,nyi-1),nx,mpi_double,precision,kdown,1,mpi_comm_world,istatus,ierr)
  call mpi_sendrecv(f(0,nyi),nx,mpi_double_precision,kdown,1,f(0,nyf+1),nx,mpi_double,precision,kup,1,mpi_comm_world,istatus,ierr)
!
  do iy=nyi,nyf
    do ix=0,nxm
      fd(ix,iy)=f(ix+1,iy)+f(ix-1,iy)+f(ix,iy+1)+f(ix,iy-1)-4.d0*f(ix,iy)
    enddo
  enddo
!
  call mpi_finalize(ierr)
  stop
end program main
解説はそのうちに