Murakami's Memorandum

Home Page > My Memo               災害科学ホワイトボード   災害科学コース   理学部

Fortran90/95入門

(工事中)

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


リンク


コンパイラー(フリーウェア)

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の値を1から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文で増分値を1ではなく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の値が2よりも小さい場合に,配列の順位と値を表示する
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 (論理式1) then
   実行文A (論理式1が真であれば実行文Aを実行する)
else if (論理式2) then
   実行文B (論理式2が真であれば実行文Bを実行する)
else
   実行文C (論理式1,2が真でない場合に実行文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

結果 4コアの場合

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

最終更新時間:2020年06月28日 19時02分49秒


HomePage > FrontPage