管理メニュー

超音波流体屋のプログラム備忘録

Fortran

最終更新：

usapfrog

- view
管理者のみ編集可
Matlabが遅いよってことで一から勉強するよ

Hello Worldまで

コンパイラだけなら、linuxは常備(なかったらyumなどで)。
[Windows]
OctaveかMingGW環境で手に入る。
Windowsならインストール後"C:\Octave\Octave3.6.4_gcc4.6.2\mingw\bin"と"C:\Octave\Octave3.6.4_gcc4.6.2\mingw\msys\1.0\bin"
または、"C:\MinGW\bin"と"C:\MinGW\msys\1.0\bin"にパスを通しとくのを忘れないこと。
(LAPACKを使う予定があるならOctaveをいれておくこと)

[Mac]
macportsやhomebrewなどで、gfortranをインストールする。 

$ vi test.f90
program main
write(*,*)'Hello World'
end
$ gfortran test.f90
$ a.exe

開発環境/IDE

まあでも、やっぱIDE欲しいよ。
(WindowsならMinGW導入後) EclipseのPhotranを導入する。
導入法に従って、インストール。
(今までのeclipseに後付しないんなら Eclipse for Parallel Application Developersを落としてきて解凍すればよし)

ついでに日本語化をしておく。
  • http://mergedoc.sourceforge.jp/の真ん中下くらいのPleiades プラグインをダウンロード・解凍
  • 中身をeclipseのフォルダにコピペ・統合する。
  • eclipse.iniの最終行に以下を追加
-javaagent:plugins/jp.sourceforge.mergedoc.pleiades/pleiades.jar
  • 初回のみ同梱のeclipse.exe -clean.cmdから起動する
  • 64bit windowsで起動時にJavaがエラーがないとエラーが出るときは、64bit javaを手動でインストールして、http://www.java.com/ja/download/manual.jsp、C:\Program Files\Java\jre7\binにパスを通すこと

  • 新規→プロジェクトでFortranぽいのを選ぶ
  • 新規→ソースでFortranソースを選ぶ
  • 上のプログラム書く
  • ビルド・実行 (合間にでるオプションは Local Fortran Application, MinGW gdb)


コード規約

流儀が本やらサイトによってそれぞれみたいだけど大体以下をベースに。
気象研究所
北海道大学 Fortran90 Tips
NAG Fortran入門
FortranTips

変数宣言

変数同士の計算式は大体C同様で計算できる。
複素数標準実装、倍精度のdの使い方が特徴的。
というよか、Fortranでは少し神経質に倍精度指定しなきゃならないらしい。
うっかり1.0とか使うと、単精度になり桁に余計なゴミが入る。 

  1. program sample
  2. integer :: a
  3. logical :: p
  4. real(8) :: y
  5. real(8), parameter :: pi = 3.14159 !定数など。(=Cのconst)
  6. complex(8) :: s
  7. character(len=4):: w
  8.  
  9. a = 3
  10. p = .true.
  11. y = 42.0d0         !dは 3e-3のeみたいに使う
  12. s = (6.0d0, 3.0d0) !これで s = 6 + 3i
  13. w = "abcd"
  14.  
  15. print *, a, z, u, w
  16. end program sample
  17.  
追記: characterはcharcter :: w*4でなく上記フォーマットのほうが了解性の点で良いらしい。
掛け算と混同してもいかんしね。

配列

足し算・内積。行数がずいぶん削れるのがありがたい。
printの書式はCとちょっと違う感じ。 
  1. program sample
  2.     implicit none
  3.  
  4.     integer, parameter :: D = 3
  5.     real(8) :: x(D) = (/1d0, 2d0,  3d0/)
  6.     real(8) :: y(D) = (/2d0, 1d0, -1d0/)
  7.     real(8) :: z(D), v(D)
  8.     real(8) :: w
  9.  
  10.     real(8) :: A(D,D)
  11.     A(1,1:3) = (/ 8.0d0, 1.0d0, 6.0d0 /)
  12.     A(2,1:3) = (/ 3.0d0, 5.0d0, 7.0d0 /)
  13.     A(3,1:3) = (/ 4.0d0, 9.0d0, 2.0d0 /)
  14.  
  15.     z(:) = x(:) + y(:)
  16.     v = matmul(A, x)
  17.     w = dot_product(x, y)
  18.  
  19.     print "(a, f8.3, f8.3, f8.3)", "z=", z
  20.     print "(a, f8.3, f8.3, f8.3)", "v=", v
  21.     print "(a, f8.3)", "w=", w
  22. end program sample
  23.  
結果 
z=   3.000   3.000   2.000
v=  28.000  34.000  28.000
w=   1.000

ループ・分岐

do i=初期,限界,増分
else if (condition) then 
  1. program sample
  2. integer :: i, j, k
  3. j=0; k=0
  4. do i=1,20
  5.   print *, "loop i=", i
  6.   k= k + 3
  7.   if ( k > 15) then
  8.     print *, "exiting loop..."
  9.     exit ! = break in C
  10.   end if
  11.   if ( k > 6) cycle ! =　continue in C
  12.   j= j + 2
  13. enddo
  14. print *, j, k
  15. end program sample
  16.  
結果 
loop i=           1
loop i=           2
loop i=           3
loop i=           4
loop i=           5
loop i=           6
exiting loop...
          4          18

関数・モジュール

追記:引数には入出力指定を示すintent(in/out/inout)を付記すること。


sampmodule.f90 
  1. module sampmodule
  2.     implicit none
  3.     real(8),parameter :: c = 340d0
  4.     real(8),parameter :: rho = 1.3d0
  5. contains
  6. integer function smd_test(a, b, x)
  7.     integer, intent(in) :: a, b, x
  8.     smd_test = a + b * x
  9. end function smd_test
  10. real(8) function smd_pres(u)
  11.     real(8), intent(in) :: u
  12.     smd_pres = rho * c * u
  13. end function smd_pres
  14. end module sampmodule
  15.  

sample.f90 
  1. program sample
  2.     use sampmodule
  3.     implicit none
  4.     integer :: d
  5.     real(8) :: p
  6.     print *,c
  7.     d = smd_test(1, 3, 5)
  8.     print *,d
  9.     p = smd_pres(2.0d0)
  10.     print *,p
  11. end program sample
  12.  


ベクトルを返すとき

戻り値がベクトルの時は、resultを指定する。例は外積関数。 
  1. function cross(a,b) result(axb)
  2.         implicit none
  3.         real(8),dimension(3) :: axb
  4.         real(8),dimension(3),intent(in) :: a
  5.         real(8),dimension(3),intent(in) :: b
  6.  
  7.         axb(1) = a(2)*b(3) - a(3)*b(2)
  8.         axb(2) = a(3)*b(1) - a(1)*b(3)
  9.         axb(3) = a(1)*b(2) - a(2)*b(1)
  10. end function cross
  11.  

テキストIO

statusのオプションが "r/w"に対応して、"old / (new, replace)"
追記をしたいときはopen文にposition='append'オプションを付ける。
エラー処理をしたければiostat=status(整数変数)をつければ、
open時読み込みエラーをstatusに格納してくれて、open時のエラーをスルーできる。 

  1. program sample
  2.     implicit none
  3.     integer :: i, j, status
  4.     integer,parameter :: fout = 11
  5.     integer,parameter :: fin = 12
  6.     real(8) :: x(20,20)
  7.     do i=1,20
  8.     do j=1,20
  9.         x(i, j) = i * 100d0 + j;
  10.     enddo
  11.     enddo
  12.     open(fout, file='mydata.dat', status='replace')
  13.     write(fout, *)  x
  14.     close(fout)
  15.     open(fin, file='mydata.dat', status='old', iostat=status)
  16.     if ( status .ne. 0 ) stop
  17.     read(fin, *)  x
  18.     close(fin)
  19. end program sample
  20.  

バイナリIO

サイズを指定して、オプションが長くなるだけ。 access, form, recl, rec 
  1. program sample
  2.     integer :: siz = 20*20*8
  3.     open(fin, file='mydata.bin', access='direct', form='unformatted', status='old', recl=siz)
  4.     read(fin, rec=1)  x
  5.     close(fin)
  6.     open(fout, file='mydata.bin', access='direct', form='unformatted', status='replace', recl=siz)
  7.     write(fout,rec=1)  x
  8.     close(fout)
  9. end program sample
  10.  

共有メモリメモリ並列(SMP)

gfortranの共有メモリ並列(-fopenmp)はコンパイラオプションの最適化(-O3)とどこか相性悪いらしく、
扱う配列サイズがやや大きい場合、segmentaion faultを吐くことがあるらしい。
この場合、最適化は(-O1)を指定するとうまくいくことがある。 
  1. program sample
  2.     implicit none
  3.     integer, parameter :: n=36
  4.     integer :: i, ene(n), res
  5.     ene(:) = 0
  6.     res = 3
  7.     !$omp parallel do private(i)
  8.     do i=1,n
  9.         ene(i) = 2*i + res
  10.     end do
  11.     !$omp end parallel do
  12.     print *, ene(:)
  13. end program sample
  14.  
バグった時はprivate指定をもう一度見なおしておくと良い。

3次逆行列/3次方程式

N>4のときは、素直にLapackとか然るべき数値計算ルーチンを使うべきだけど、
多用する3x3くらいは、(必要なライブラリかき集めるのめんどいし、四角い車輪の再生産にあるけど)
自分で書いていい気がする。 
  1. function inv3(A) result(Y)
  2.         implicit none
  3.         real(8), intent(in) :: A(3,3)
  4.         real(8) :: Y(3,3)
  5.         real(8) :: det
  6.         integer :: i,j
  7.         det = 0
  8.         do i=1,3
  9.             det = det + A(i,1)*cofact3(A,i,1)
  10.         end do
  11.         do i=1,3
  12.             do j=1,3
  13.                 Y(i,j) = cofact3(A,j,i) / det
  14.             end do
  15.         end do
  16. end function inv3
  17. real(8) function cofact3(A,i,j)
  18.         implicit none
  19.         real(8), intent(in) :: A(3,3)
  20.         integer, intent(in) :: i, j
  21.         integer :: i2, i3, j2, j3
  22.         i2 = mod(i, 3)+1; i3 = mod(i+1, 3)+1;
  23.         j2 = mod(j, 3)+1; j3 = mod(j+1, 3)+1;
  24.         cofact3 = A(i2,j2)*A(i3,j3)-A(i2,j3)*A(i3,j2)
  25. end function cofact3
  26.  
  27.     !wikipediaの三次方程式の解の公式・実部しか取ってこないので注意
  28. subroutine solve3(x, a)
  29.         implicit none
  30.         real(8), intent(out) :: x(NDIM)
  31.         real(8), intent(in) ::  a(NDIM)
  32.         integer :: k
  33.         complex(8) :: j, w, y(NDIM), p, q, R, S
  34.         j = (0, 1)
  35.         w = 0.5d0*(-1d0+j*sqrt(3d0))
  36.         p = a(2)-a(1)**2/3d0
  37.         q = a(3)-a(2)*a(1)/3d0+2d0/27d0*a(1)**3
  38.         R=(-q/2d0 + sqrt((q/2d0)**2 + (p/3d0)**3))**(1d0/3d0)
  39.         S=(-q/2d0 - sqrt((q/2d0)**2 + (p/3d0)**3))**(1d0/3d0)
  40.         do k=1,3
  41.             y(k) = w**(k-1) * R + w**(4-k) * S
  42.             x(k) = real(y(k),8) - a(3)/3d0
  43.         end do
  44. end subroutine solve3
  45.  


「Fortran」をウィキ内検索
人気記事ランキング
  1. バターワースフィルタ定数
  2. Lapack
  3. CLapack
  4. NLopt
  5. Linux計算サーバー構築メモ
もっと見る
最近更新されたページ
もっと見る
目安箱バナー