2010年1月9日土曜日

[Fortran] if文とstop文

以下のプログラムは,1から100までの整数の内井,偶数と奇数の和をそれぞれ計算するプログラムです.

program loop_odd_even
  implicit none
  integer i, sum0, sum1
  sum0 = 0
  sum1 = 0
  do i = 1, 100
    if (mod(i, 2) == 0) then
      sum0 = sum0 + i
    else if (mod(i, 2) == 1) then
      sum1 = sum1 + i
    else
      stop 'sommething is wrong !!'
    endif
  enddo
  write(*, *) 'sum0, sum1, sum =', sum0, sum1, sum0 + sum1

endprogram loop_odd_even

doループの内部では,制御変数 i が偶数か奇数かを判別するために,組み込み関数である余り関数(mod(a, p))をif文に利用して,条件判定を行なっています.
このプログラムをコンパイル,実行すると以下のようになります.
$ gfortran loop_odd_even.f90
$ ./a.out
 sum0, sum1, sum =        2550        2500        5050


if文で最も多く用いられるのは,数値を比較する論理式で,以下のような関係演算子を用いて記述されます.

if (a > b)    abより大きければ真
if (a >= b)   abより大きいか等しいとき真
if (a <= b)   abより小さければ真
if (a == b)   abより小さいか等しいとき真

if (a /= b)   abと等しくなければ真

以下のプログラムは,上記のプログラムではプログラム中で指定していた合計を計算する範囲をキーボード入力(標準入力)で定めるように書き換えたものです.

program loop_inf
  implicit none
  integer i, n, sum

  do
    write(*,*) ' input n (if n <= 0, stop) :'
    read(*,*) n
    if (n <= 0) stop ' good bye ...'
    sum = 0
    do i = 1, n
      sum = sum + i
    end do 
    write(*,*) 'sum = ', sum
  enddo

end program loop_inf

$ gfortran loop_inf.f90
$ ./a.out
  input n (if n <= 0, stop) :
10
 sum =           55
  input n (if n <= 0, stop) :
100
 sum =         5050
  input n (if n <= 0, stop) :
0
STOP  good bye ...

上記のプログラムの中では,read文が使われています.
read( )の括弧内には,write( )の括弧と同様に,ファイル番号と書式をカンマで区切って記述します.
read文では,ファイル番号は入力元を指定し,これを"*"とすると,標準入力(処理系で定められた入力元で,通常はキーボード)からの入力となります.
また,read文の書式を"*"とすると並び入力となり,read文の後に並ぶ変数の型に合った書式で読み込みが行われます.
上記のプログラムでは,整数型変数nと同じ型のデータとして読み込みを行います.このため,数字以外の英字や記号を入力するとエラーとなります.また,実数(ex 1.23)を入力すると小数点以下が切り捨てられた整数とみなされます.例えば,0.01と入力すると,切り捨てにより n = 0 となり,上記のプログラムでは処理が終了となります.

2010年1月8日金曜日

[Fortran] DOループの処理速度

DOループにおいて,ソースコードの書き方によってアクセスに要する時間が異なります.

整数型2次元配列のという"nteger(300090000)"各要素に,行番号(j)と列番号(k)の和を代入する例を用いて,計算処理速度の違いを比較してみます.

以下のプログラムの中で,組込み関数"date_and_time"は,計算処理時刻を知る内部服プログラムです.

PROGRAM TimeRace
IMPLICIT NONE
  CHARACTER:: date*8, time*10, zone*5   !Character data
  INTEGER::   nval(8)                   !Integer type
  INTEGER::   j, k
  INTEGER::   nteger(3000, 90000)
! Slow access example
  CALL date_and_time(date, time, zone, nval)
  WRITE(*,*)  time
    DO J=1, 3000
      DO k=1, 9000
        nteger(J, k)=j+k
      ENDDO
    ENDDO
  CALL date_and_time(date, time, zone, nval)
    WRITE(*,*) time
! Early access example
  DO k=1, 90000
    DO j=1, 3000
      nteger(J, k)=j+k
    ENDDO
  ENDDO
  CALL date_and_time(date, time, zone, nval)
    WRITE(*,*)  time
END PROGRAM TimeRace

このプログラムをコンパイルして実行すると,以下のようになります.

$ gfortran Timerace.f90
$ ./a.out
 095240.164
 095250.816
 095251.825
$ ./a.out
 095255.929
 095305.783
 095306.793

比較のために2回実行していますが,1回目の実行では,まず最初のDOループをスタートした時間
 095240164(9時52分40.164秒)
が表示され,続いて最初のDOループが終了した時間
 095250816(9時52分50.816秒)
が表示されます(最初のDOループでは,処理に約10秒かかっています).
続いて,2番目のDOループが終了した時間
 0952251.825(9時52分51.825秒)
が表示されます(2番目のDOループでは,処理に約1秒かかっています).
確認のための2回目の実行でも,1回目は約10秒,2回目は約1秒の処理時間です.

一般的に,外側ループの制御変数(配列の右側添字)をある値に固定したままで,内側ループにある制御変数(配列の左側添字)を変化させるようなDOループが効率的とされており,このことを確認することができました.

*配列について*
Fortran90/95における配列は,数学における行列と同じような並びとなります.

のような行列は,

となります.
すなわち,行列Aにおける1行2列成分(a_12)は,配列Aの1行2列成分(a(1, 2))と同じです.

2010年1月7日木曜日

[Fortran] 処理時間を計測する cpu_time

処理時刻を知る組み込み関数として,組み込み関数"cpu_time"を用いることができます.
"cpu_time"は,Fortran90では一部の処理系のみで使用可能で,Fortran95では標準として組み込まれています.
この"cpu_time"を用いて処理作業を挟む2箇所で時間カウンタ(引数)を取り,両カウンタの差から処理時間を秒単位で評価することができます."cpu_time"を用いるには,

 call cpu_time(time)

で呼び出します.ここに,"time"は実数型の引数で,"intent(out)"の属性を持ちます.
意味のある時間を返せない場合は,負値が帰ってきます.プログラム例を以下に,示します.

program cpu_time
  implicit none
  real(8):: t1, t2, t_req
  call cpu_time(t1)
!
! arithmetic processing
!
  call cpu_time(t2)
    t_req = t2 - t1 !time required
    write(*,*) 'Time required in sec = ', t_req
end program cpu_time

このプログラムを"cpu_time.f90"という名前でコンパイル,実行すると以下のようになります.
$ gfortran cpu_time.f90
$ ./a.out
 Time required in sec =    2.9999999999995308E-006
このプログラムでは,演算処理部分が空なので,"cpu_time(t1)"と"cpu_time(t2)"の間の時間は"2.9999999999995308E-006 (s)"となっています.
演算処理部分に以下のように記述したプログラムを"cpu_time2.f90"という名前で保存します.

program cpu_time2
  implicit none
  real(8):: t1, t2, t_req
  integer i, a, b, c
  a = 2
  b = 1

  call cpu_time(t1)

    write(*,*) 'n = ', b, b
    write(*,*) 'n = ', a, a

    do i =3, 100000
      c = a + b
      write(*,*) 'n = ', i, c
      b = a
      a = c
    end do

  call cpu_time(t2)

    t_req = t2 - t1 !time required
    write(*,*) 'Time required in sec = ', t_req
end program cpu_time2

コンパイル実行すると以下のようになります.
$ gfortran cpu_time2.f90
$ ./a.out
 n =            1           1
 n =            2           2
………………………………………………………………………………………………………………………
abridgement
………………………………………………………………………………………………………………………
 n =        99999   873876091
 n =       100000   679394397
 Time required in sec =   0.31460899999999997     
この実行例を見ると,"cpu_time(t1)"と"cpu_time(t2)"の間の時間は"0.31460899999999997(s)"となっており,先の結果に比べて時間がかかっていることがわかります.

2010年1月6日水曜日

[Fortran] 配列を使った内積を計算するプログラム

ファイルから2次元平面内の位置ベクトルu, vの成分(u1, u2), (v1, v2)を読み取り,2つのベクトルの内積を計算するプログラム例を以下に示します.

program DP
  implicit none
  integer i
  real(8) u(2), v(2), dotp

  open(10,file = 'input1.d')
  read(10,*) (u(i), i = 1, 2) !改行の有無に関係なく配列要素の値を読み込む場合
  close(10)

  open(11,file = 'input2.d')
  read(11,*) (v(i), i = 1, 2) !改行の有無に関係なく配列要素の値を読み込む場合
  close(11)

   write(*,*) 'u_i = ', (u(i), i = 1, 2)
   write(*,*) 'v_i = ', (v(i), i = 1, 2)

   dotp = 0.0d0
   do i = 1, 2
     dotp = dotp + u(i) * v(i)
   enddo

   write(*,*) 'dot product = ', dotp

   open(20,file='output.d')
    write(20,'(2e20.4)') , u(1), u(2)
    write(20,'(2e20.4)') , v(1), v(2)
    write(20,'(e20.4)'), dotp
   close(20)

end program DP

このプログラムをターミナルからコンパイル,実行すると以下のようになります.
$ gfortran DP.f90
$ ./a.out
 u_i =    1.2000000000000000        3.3999999999999999     
 v_i =    4.0999999999999996        2.6000000000000001     
 dot product =    13.759999999999998     
ターミナルでの実行結果とは別に,'output.d'というファイルが作られ,その中身は,以下のようになります.
          0.1200E+01          0.3400E+01
          0.4100E+01          0.2600E+01
          0.1376E+02

--プログラムの簡単な説明--
プログラム中では,配列を利用して,ベクトルの成分を以下のようにあらわしています.
  u1, u2 → u(1), u(2)
  v1, v2 → v(1), v(2)
矢印(→)の右側にあるのは,配列要素(array element)です.

ファイルから成分の値を読み取る際には,改行から1つずつ配列要素の値を読み込む場合と,改行の有無に関係なく配列要素の値を読み込む場合があります(上記のプログラムは後者の例).
改行から一つずつ読み込む場合は,ベクトルuについては

  read(10,*) (u(i), i = 12!改行の有無に関係なく配列要素の値を読み込む場合

の部分を,以下のように書き換えます.

  do i = 1, 2        !改行から1つずつ配列要素の値を読み込む場合
    read(10,*) u(i)
  enddo

ベクトルvについても同様で,

  read(11,*) (v(i), i = 1, 2) !改行の有無に関係なく配列要素の値を読み込む場合

の部分を

  do i = 1, 2        !改行から1つずつ配列要素の値を読み込む場合
    read(11,*) v(i)
  enddo

に書き換えます.

read文で書式を指定すると,指定の書式に合わない入力データは読み込めなくなるので,プログラム'DP'のように,並び入力とするのが一般的です.
並び入力では,データの形式に依存せずに変数の型に合わせて入力が行われます.
同一行に書かれたデータ間の区切りとしては,1つ以上のスペース,あるいはタブを入れておけば区切りと見なされます.

ディスプレイ上にも出力した上で,ファイルにも書き出しを行うようにしており,ファイルに出力する際には,'(2e20.4)',  '(e20.4)'として,出力書式を指定しています.
'2e20.4'は,全体の文字数を20個確保して,指数形式で小数点以下4桁までが表示されるようにする書式で,2つの出力が改行なし(並び)で出力されます.
また,'e20.4'は,文字数20を確保して,指数形式で小数点以下4桁までを表示する書式です.

2010年1月5日火曜日

[Fortran] 入力の書式関数

入力の際に,出力と同様に書式を指定することができます.
ただし,書式を指定すると,その通りに入力(キーボード,もしくはファイル)の記述を行わないと正しく読み込みが行われません.
通常は,データを1つ以上のタブやスペースで区切って入力し,プログラム中では'READ(ファイル番号,*)'として並び入力により読み込みを行うことが一般的です.

書式を指定する場合には,'*'の部分に,書式付き出力の場合と同様の書式を記述します(書式を指定する場合は,連続して記述された文字や数字を指定した桁数で切り分けて読み込む場合が多い).

また,文字列を読み込む際に,書式を使うか否かによって違いが生じます.
以下のような内容の入力ファイル'data.d'があったとします.

abcde fghij
ABCDE FGHIJ

このファイルを,次のようなプログラムで読み込むとします.

PROGRAM ReadCharacter
  IMPLICIT NONE
  CHARACTER(10) char1, char2
  OPEN(10, file = 'data.d')
  READ(10, *) char1
  READ(10, '(a)') char2
  WRITE(*,*) TRIM(char1)
  WRITE(*,*) TRIM(char2)

END PROGRAM ReadCharacter

このプログラムをコンパイル,実行すると以下のようになります.

$ gfortran ReadCharacter.f90
$ ./a.out
 abcde
 ABCDE FGHI

出力結果からわかるように,'READ(10, *) char1'とすると,行中の空白データが区切りとみなされ,'abcde'の5文字が読み取られます.一方で,'READ(10'(a)') char2'と書式を指定すると,空白も含めた文字数分の文字が読み込まれます.

2010年1月4日月曜日

[Fortran] 出力の書式関数

入出力の際に数値や文字を整形して体裁を整えるには,書式を指定します.
以下に,書式付き出力の例を示すためのプログラム例を示します.

PROGRAM format
  CHARACTER(30):: CHAR1 = 'Fortran is'
  CHARACTER(20):: CHAR2 = ' FOrmula TRANslation.'
  WRITE(*,'(a7)') CHAR1
  WRITE(*,'(a,a)') CHAR1, CHAR2
  WRITE(*,'(a,a)') TRIM(CHAR1), TRIM(CHAR2)

END PROGRAM format

このプログラムをコンパイル,実行すると以下のようになります.

$ gfortran format.f90
$ ./a.out
Fortran
Fortran is                     FOrmula TRANslation.
Fortran is FOrmula TRANslation.

上記のように,write文の括弧内のカンマの右側に')'という区切りを設けて,その括弧内に書式を記述します.
文字型変数の場合には,'(a*)'という書式を用い,'*'の部分には整数(変数名ではなく数字)を書きます.そうすると,文字列の先頭からその整数の個数分の文字が出力されます.なお,整数を省略すると,変数の全ての文字が出力されます.
上記の例では,'WRITE(*,'(a7)') CHAR1'として文字列'CHAR1'の先頭の7文字分を出力しています.
また,5行目のようにWRITE(*,'(a,a)') CHAR1, CHAR2とすると,2つの文字型変数の後続の空白部分も含む全ての文字列が出力されます(3行目で'CHARACTER(20):: CHAR2 = ' FOrmula TRANslation.''として,'FOrmula'の前に半角の空白があることに注意).

組込み関数を使用すると,後続の空白を取り除いた文字列が出力されます .

また,以前の投稿にもあるように,キーボード入力を促す際には,以下の出力を行います.
      WRITE(*,fmt='(a)', advance='no''Input n: '
書式の'(a)'出力対象はWRITE文の後に続く文字定数'Input n: 'です.
また,'advance='no''は,改行を抑制する指定です.コンパイラによっては,
      WRITE(*,'(a)', advance='no''Input n: '
      WRITE(*,'(a\)''Input n: '
としても良いようですが,gfortranではコンパイルエラーとなります.ちなみに,デフォルトでは,'advance='yes''という設定になっており'WRITE'文は1回実行される毎に改行を行う設定となっています.,'advance='no''は,改行を抑制する指定です.

次に,整数型,および実数型の変数に対する書式付き出力プログラム例を示します.

PROGRAM format
  IMPLICIT NONE
  INTEGER:: i, ivar(4) = (/ (i, i = 1001, 1004) /)
  REAL(8):: xvar(4) = (/ 1.234d1, 0.0d0, -1.234d-2, -9.999d-1 /)

  DO i = 1, 4
    WRITE(*,'(i10, f10.2, e12.4)') ivar(i), xvar(i), xvar(i)
  ENDDO

END PROGRAM format

このプログラムをコンパイル,実行すると以下のようになります.
$ gfortran format.f90
$ ./a.out
      1001     12.34  0.1234E+02
      1002      0.00  0.0000E+00
      1003     -0.01 -0.1234E-01
      1004     -1.00 -0.9999E+00
このプログラムでは,整数型変数ivar(i)の出力書式に関して,'i10'という書式が指定されています.これは,全体の文字数を10個分確保して,その中に右詰めで整数値を書き出すという書式です.もし,全体の文字数が足りないとオーバーフローを起こして,文字数分の*が出力されて,値は確認できません(以下の例は.'i2'とした例).
**     12.34  0.1234E+02
**      0.00  0.0000E+00
**     -0.01 -0.1234E-01
**     -1.00 -0.9999E+00
続いて,''は,全体の文字数を10個確保して,その中で小数点以下2桁の数値を右詰めで書く書式です.
全体の文字数をw,小数点以下の数値の個数をdとすると,小数点と1の位の数値,また,マイナス符号が付く文字数を考慮すると,w≧d+3であることが最低限必要になります.
例えば,-10を'f5.2'という書式で出力しようとすると,-10.00と文字数が6個必要になるので,オーバーフローを起こしてしまいます.

数値計算では,上記のような'f5.2'という形式よりも,'e12.4'のような指数形式の書式(出力)が使われる場合が多くなります.'e12.4'は,全体の文字数を12個確保して,10の指数形式(12.34 を 0.1234x10^2と表す形式)で,小数点以下を4桁表示する書式です.
全体の文字数と小数点以下の数値の個数をそれぞれw,およびdと表すとき,w≧d+7であることが必要となります(指数の符号,数値2桁分,E,先頭の符号,1の位の数値,小数点の文字数が必要であるため).

同様な書式を繰り返し利用するには,以下のような書式指定があります.
    WRITE(*,'(4e12.4)') xvar(1:4)
この書式では,同じ行に'e12.4'の書式で4つの出力が行われます.

2010年1月3日日曜日

[Fortran] 無限ループ

Fortran 90/95において,キーボード入力した値までの自然数の和を求めるプログラム例を以下に示します.

PROGRAM loop_input
  IMPLICIT NONE
  INTEGER sum, n, i
    DO
      WRITE(*,*) 'Input n (if n <= 0, Stop) : '
      READ(*,*) n
      if (n <= 0) STOP 'end'
      sum = 0
      DO i = 1, n
        sum = sum + i
      ENDDO
      WRITE(*,*) 'Sum = ', sum
    ENDDO
END PROGRAM loop_input

上記のプログラムでは,入力と結果の表示が繰り返し行われるように,制御変数なしのDOループとしています.
また,プログラムにおいては,2つのDOループが使われていて,一方が他のループに含まれる形になっており,これを2重ループといいます.外側のDOループは,制御変数なしのDOループなので,DO - ENDDOの間にループから抜け出す仕組みを作っておかないと,反復演算が無限に繰り返されることになります.上記のプログラムでは,入力値が0以下であればSTOP文が実行されて終了するようになっています.

このプログラムをコンパイル,実行すると以下のようになります.

$ gfortran loop_input.f90
$ ./a.out
 Input n (if n <= 0, Stop) : 
50
 Sum =         1275
 Input n (if n <= 0, Stop) : 
100
 Sum =         5050
 Input n (if n <= 0, Stop) : 
-3
STOP end

上記の例では,キーボードからの入力を促す表示の後に改行されていますが,この改行を行わない場合は,WRITE文を以下のように書き換えます.

      WRITE(*,fmt='(a)', advance='no') 'Input n (if n <= 0, Stop) : '

このように書き換えてコンパイル,実行すると以下のようになります.
$ gfortran loop_input.f90
$ ./a.out
Input n (if n <= 0, Stop) : 30
 Sum =          465
Input n (if n <= 0, Stop) : 88
 Sum =         3916
Input n (if n <= 0, Stop) : 0
STOP end