Fortran90/95入門

（工事中）

　学部の3年生の必修の授業でならった言語がFortranで，それ以来お世話になってます。最近は大抵のことはMatlab/Octave/Rで済ませますが，計算速度の必要な計算はやはりFortranです。というわけで，説明用の資料として。

• リンク
• コンパイラー（フリーウェア）
• 超入門編
  • 最初の一歩
  • とりあえずひな形
  • 繰り返し(do loop)
  • 条件分岐（if 文）
• 中級編
  • 配列変数の演算(Fortran90/95の魅力の一つ)
  • 副プログラム（サブルーチン）
  • モジュール
• 上級編　
  • gfortran + Open MP 並列化
    • サンプルプログラム
    • コンパイル　　
    • 結果　４コアの場合
  • ifortran + Open MP
    • コンパイル
  • FFTW3を使う
• いろいろメモ

リンク

• データ解析のためのFortran90/95入門 IIAE（一般財団法人大気環境総合センター） 「データ解析のためのFortran90/95入門　愛媛大学農学部（気象系）」を収録
• Fortran標準コーディングルール　気象研究所 気象関係者間でのソースコードの流通を目的としたコーディングルール
• Fortran演習（地球惑星物理学演習）　東京大学
• 2. Fortran90/95入門 神戸大学　陰山先生
  • 2.1 Fortran90/95入門前半
  • 2.2 Fortran90/95入門後半
• 出版年の新しい書籍紹介
  • 牛島　省，『数値計算のためのFortran90/95プログラミング（第2版）』，森北出版(2020/1/28)
    • 追加・訂正 　2019年5月時点でのインストールについての推奨

コンパイラー（フリーウェア）

• GNU Fortran　gfortran
• The G95 Project　g95

Intel Fortran (ifort) 非商用版

超入門編

 最初の一歩

プログラム・ファイルの名前（拡張子）の付け方

• fortran90/95規格のソース　sample.f90
• fortran77規格のソース　sample.f

プログラムのコンパイル

• gfortran sample.f90 -o sample
• g95 sample.f90 -o sample
• ifc sample.f90 -o sample

 とりあえずひな形

program sample　　! プログラム文（プログラム名を宣言）
  implicit none　! 使用する全ての変数の型宣言をするという宣言
  integer i, j　　! 整数を宣言
  real(8) x       ! 実数（倍精度）を宣言
  ...
  ...
end program sample　　! end文

• プラグラム文は省略可能であるが，end文は必須。プログラム文を省略した場合には単にendのみ。
• real(8)は，実数を倍精度（約15桁）の変数にするという宣言。例えば，Excel等では標準でこの程度の精度で演算をしているので，メモリを節約したい時以外にはあえて基本形real（約7桁）を使わなくても良い。

 繰り返し(do loop)

• do〜enddo の中の処理をiの値を１から10まで変えながら繰り返す。

program sum
  implicit none
  integer i
  real(8) sum, x(10)　　! 実数変数sumと実数配列変数x（要素数が10）を宣言
  do i=1,10
     x(i)=0.5*i　　! 配列変数に値を設定する　x(1)=0.5, x(2)=1.0, ..., x(10)=5.0
  enddo
  sum=0
  do i=1,10
     sum=sum+x(i)　　! 合計x(1)+x(2)+...+x(10)を計算する
  enddo
  write(*,*) 'sum= ',sum
end program sum

• fortran77 では行番号を使い繰り返しの範囲を示すが，fortran90/95では行番号を使用しない。

      do 10 i=1,10
10    x(i) = 0.5*i

• do~enddo文で増分値を１ではなく2とする場合 (i=1,3,5,7,9となる)

do i=1,2,10
  x(i)=0.5*i
enddo

• 無限ループ　（抜けるための対策が必要：　if, exit, stop, goto 文）

do
  x(i)=0.5*i
enddo

 条件分岐（if 文）

• 配列変数mの値が２よりも小さい場合に，配列の順位と値を表示する

program sample1
implicit none
integer :: m(5)=(/-1, 5, 3, -4, 2/) ,i   ! m(1)=-1, ..., m(5)=2 を初期値として設定
  do i=1,5
    if (m(i)<2) then
      write(*,*) ' m<2', i, m(i)
    endif
  enddo
end program sample1

• if文の基本的な構造

if (論理式１) then
   実行文A　（論理式１が真であれば実行文Aを実行する）
else if (論理式２) then
   実行文B　（論理式２が真であれば実行文Bを実行する）
else
   実行文C　（論理式１，２が真でない場合に実行文Cを実行する）
endif

• 論理式の構造

if (a==b)   a と b が等しければ真
if (a>b)     a が b より大きければ真
if (a>=b)   a が b 以上であれば真
if (a<b)     a が b よりも小さければ真
if (a<=b)   a が b 以下であれば真
if (a/=b)     a と b が等しくなければ真

中級編

 配列変数の演算(Fortran90/95の魅力の一つ)

!  do i=1,10
!      x(i)= y(i)*5+1
!  enddo

   x(1:10)=y(1:10)*5+1

 副プログラム（サブルーチン）

program main
    implicit none
    integer i,num
    real(8) sum
    ...
    call sub1(num,sum)
    ...
end program main
subroutine sub1(num,sum)
     implicit none
     integer i,num
     real(8) sum
     ...
end subroutine sub1
  

 モジュール

fortran90以降で出てきた新しい規格
（工事中）

上級編　

  gfortran + Open MP 並列化

gfortranはデフォルトで並列化に対応していました！

Fortran プログラマのためのOpenMP 入門　京都大学学術メディアセンター
OpenMP入門

サンプルプログラム

program testOMP2
  !$ use omp_lib
  implicit none
  print *, 'START'
!$omp parallel
  print *, 'Hello! N =', omp_get_num_threads(), ' and I am ', omp_get_thread_num()
!$omp end parallel
  print *, 'END'
end program testOMP2

コンパイル　　

% gfortran -fopenmp testomp2.f90

結果　４コアの場合

start
Hello! N=        4 and I am        0
Hello! N=        4 and I am        1
Hello! N=        4 and I am        2
Hello! N=        4 and I am        3
END

 ifortran + Open MP

インテルFortranコンパイラーOpenMP活用ガイド

コンパイル

ifort -openmp -openmp-report2 test_omp.f90

 FFTW3を使う

1. FFTW3をダウンロード
2. コンパイル
   1. ./configure
   2. make
   3. make install
3. fftw3を使用するソースをコンパイルする

% gfortran tfftw3.f90 -I/usr/local/include -L/usr/local/lib -lfftw3

ubuntuのリポジトリのlibfftw3やlibfftw3-devを使用する場合
% gfortran tfftw3.f90 -I/usr/include -lfftw3 -o tfftw3
 (ライブラリやインクルードファイルのディレクトリが異なる)

• 計算例

sample program

program fftw3_test
implicit none
! array size
integer, parameter :: N=16
! input data
real(8), dimension(N) ::in,in2
! output data 
complex*16, dimension(N/2+1) :: out
! plan
integer(8) ::plan

integer    ::i,M
include 'fftw3.f'

!=========================
! set data
!========================= 
in(1)=5.0
in(2)=32.0
in(3)=38.0
in(4)=-33.0
in(5)=-19.0
in(6)= -10
in(7)= 1.0
in(8)= -8.0
in(9)= -20.0
in(10)= 10.0
in(11)= -1.0
in(12)= 4.0
in(13)= 11.0
in(14)= -1.0
in(15)= -7.0
in(16)= -2.0
in2(1:N)=in(1:N)
!=====================
call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(plan,in,out)
call dfftw_destroy_plan(plan)
!=====================
   M=N/2+1
do i=1,M
   print *,real(out(i)),aimag(out(i)) 
end do
! 
!=====================
call dfftw_plan_dft_c2r_1d(plan,N,out,in,FFTW_ESTIMATE)
call dfftw_execute_dft_c2r(plan,out,in)
call dfftw_destroy_plan(plan)
!=====================
do i=1,N
   print *, in(i)/dble(N),in2(i)
end do
end program fftw3_test

入力が実数の場合のFFT演算では、出力のメモリサイズはN/2+1とする。半分は出力されない。演算結果は通常のFFT演算（Octave等）と同じ。また、逆変換でもメモリサイズの関係はかわらない。

演算結果
　murakami@murakami-VB:~/work$ ./tfftw3

   0.0000000000000000        0.0000000000000000     
   62.073801020876488        33.141504861910015     
   43.911688245431428       -67.041630560342611     
   39.663452799755568       -95.619529124742456     
  -54.000000000000000       -70.000000000000000     
  -33.504073233321513        30.848508306793022     
  -57.911688245431428        18.958369439657382     
   31.766819412689461        39.609542293445493     
   16.000000000000000        0.0000000000000000     <----- ここまでがフリエ変換の出力
   5.0000000000000000        5.0000000000000000     <----- 逆変換された値は規格化する（1/N)
   32.000000000000000        32.000000000000000     
   38.000000000000000        38.000000000000000     
  -33.000000000000000       -33.000000000000000     
  -19.000000000000000       -19.000000000000000     
  -10.000000000000002       -10.000000000000000     
   1.0000000000000033        1.0000000000000000     
  -7.9999999999999982       -8.0000000000000000     
  -20.000000000000000       -20.000000000000000     
   10.000000000000002        10.000000000000000     
  -1.0000000000000040       -1.0000000000000000     
   3.9999999999999973        4.0000000000000000     
   11.000000000000000        11.000000000000000     
  -1.0000000000000013       -1.0000000000000000     
  -7.0000000000000018       -7.0000000000000000     
  -1.9999999999999984       -2.0000000000000000     

いろいろメモ

ubuntu10.04LTS/Intel fortran 11.0

• アポストロフィーを含む文字列の出力(f77)

character*80 text1
text1="set dir=''0 0.01 0''"
L1=len_trim(text1)
write(10,'(a<L1>)') text1

ダブルコーテーション(")中のアポストロフィー(')は二重にする。

set dir='0 0.01 0'

f90でのメモ

program gtest
implicit none
character ::text1*80,text2*80,ciday*80
integer iday,L1
iday=30
text1="set d='0 0.1 0.2'"
text2='opaq url'
write(*,'(a)') text1
write(*,'(a)') text2
write(ciday,*) iday
ciday=adjustl(ciday)
write(*,'(a)') trim(text2)//trim(ciday)//trim(text1)
L1=len_trim(text2)
write(*,'(a)') text2(1:L1)
end program gtest