Fortran90入門


FORTRAN77を学習してプログラミングできるようになった気がしますので、Fortran90の方へ移っていきたいと思います。

Fortran90 ビルド、継続行

Fortran90をビルドするには、拡張子を「.f」から「.f90」に変更すれば良いようです。しかも77の仕様で作成したソースコードは、Fortran90で使えるようです。
では、試してみます。

takk@deb9:~$ gfortran cut8.f -fdollar-ok

77でビルド後、90用にファイル名を変更して、ビルドします。

takk@deb9:~$ mv cut8.f cut8.f90
takk@deb9:~$ gfortran cut8.f90 -fdollar-ok
cut8.f90:12:56:

             OPEN(UNIT=10,FILE=ARGS(1:INDEX(ARGS,' ')-1),
                                                        1
Error: Syntax error in OPEN statement at (1)
cut8.f90:13:12:

      1      STATUS='OLD',
            1
Error: Unclassifiable statement at (1)
cut8.f90:14:12:

      2      ACCESS='SEQUENTIAL')
            1
Error: Unclassifiable statement at (1)
takk@deb9:~$

あれ、ビルド通りません。行の連結が上手くいっていないのでしょうか。
3行に渡ってOPEN関数を書いていましたが、1行にまとめてみます。

            OPEN(UNIT=10,FILE=ARGS(1:INDEX(ARGS,' ')-1),
     1      STATUS='OLD',
     2      ACCESS='SEQUENTIAL')

これが、こうなりました。

            OPEN(UNIT=10,FILE=ARGS(1:INDEX(ARGS,' ')-1),STATUS='OLD',ACCESS='SEQUENTIAL')
 

ビルドできました。

takk@deb9:~$ gfortran cut8.f90 -fdollar-ok
takk@deb9:~$

無理やり一行にしてしまいましたが、Fortran90で継続行の書き方が変わっているようですので、それに対応してみます。
行が継続する場合は、&を行末につけて、その次の行の開始にも&が必要になります。

            OPEN(UNIT=10,FILE=ARGS(1:INDEX(ARGS,' ')-1), &
     &      STATUS='OLD', &
     &      ACCESS='SEQUENTIAL')

ビルドして確認。

takk@deb9:~$ gfortran cut8.f90 -fdollar-ok
takk@deb9:~$

通りました。
&の位置はそんなに厳密でもないようです。

            OPEN(UNIT=10,FILE=ARGS(1:INDEX(ARGS,' ')-1), &
           &STATUS='OLD', &
           &ACCESS='SEQUENTIAL')
takk@deb9:~$ gfortran cut8.f90 -fdollar-ok
takk@deb9:~$

ずっと気になってますが、Fortran90になっても、-fdollar-okって必要なのでしょうか。はずしてビルドしてみます。

takk@deb9:~$ gfortran cut8.f90
cut8.f90:54:16:

  3300   FORMAT(A$)
                1
Fatal Error: Invalid character ‘$’ at (1). Use ‘-fdollar-ok’ to allow it as an extension
compilation terminated.
takk@deb9:~$

エラーになったので必要そうです。

Fortran90 do

今まで繰り返し制御に使っていたDO文は、こんな形でした。

takk@deb9:~$ cat do.f
      PROGRAM MAIN
        DO 1000 I=1,10
          WRITE(*,*) 'HELLO'
 1000   CONTINUE
        STOP
      END
takk@deb9:~$ gfortran do.f
takk@deb9:~$ ./a.out
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
takk@deb9:~$

手元の古い教科書がCONTINUEを使う構文でしたので気にしてませんでしたが、どうやら古い書き方のようで、繰り返しは、END DOを使うのが正解のようです。

takk@deb9:~$ cat do.f90
      PROGRAM MAIN
        DO I=1,10
          WRITE(*,*) 'HELLO'
        END DO
        STOP
      END
takk@deb9:~$ gfortran do.f90
takk@deb9:~$ ./a.out
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
takk@deb9:~$

END DOでも行番号を指定して、入れ子のDOに対して一つのEND DOにできるのでしょうか。
入れ子にする前に、行番号が使えるかどうかやってみます。

takk@deb9:~$ cat do.f90
      PROGRAM MAIN
        DO 10000 I=1,10
          WRITE(*,*) 'HELLO'
10000   END DO
        STOP
      END
takk@deb9:~$ gfortran do.f90
takk@deb9:~$ ./a.out
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
takk@deb9:~$

使えるようです。
次に入れ子にします。ただし、DO 2つに対して、END DOは1つです。

takk@deb9:~$ cat do.f90
      PROGRAM MAIN
        DO 10000 I=1,2
         DO 10000 J=1,3
          WRITE(*,*) 'HELLO'
10000   END DO
        STOP
      END
takk@deb9:~$ gfortran do.f90
takk@deb9:~$ ./a.out
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
 HELLO
takk@deb9:~$

許容されてるんですね。
では、行番号を外すとどうなるのでしょう。

takk@deb9:~$ cat do.f90
      PROGRAM MAIN
        DO I=1,2
         DO J=1,3
          WRITE(*,*) 'HELLO'
        END DO
        STOP
      END
takk@deb9:~$ gfortran do.f90
do.f90:7:9:

       END
         1
Error: END DO statement expected at (1)
f951: Error: Unexpected end of file in ‘do.f90’
takk@deb9:~$

ビルドエラーになりました。
使い方は分かりました。CONTINUEを使う書き方は古いらしいので、これからはEND DOを使っていきます。

Fortran90 do exit

Fortran、次は、do構文のexitを使ってみます。

まずはexitなしの繰り返し文です。

takk@deb9:~$ cat do-only.f90
      PROGRAM MAIN
        DO I=1,2
          DO J=1,3
            WRITE(*,*) I,J
          END DO
        END DO
        STOP
      END
takk@deb9:~$ gfortran do-only.f90
takk@deb9:~$ ./a.out
           1           1
           1           2
           1           3
           2           1
           2           2
           2           3
takk@deb9:~$

EXITを追加します。
Jが2になったらEXITする。

takk@deb9:~$ cat do-exit.f90
      PROGRAM MAIN
        DO I=1,2
          DO J=1,3
            IF(J.EQ.2)EXIT
            WRITE(*,*) I,J
          END DO
        END DO
        STOP
      END
takk@deb9:~$ gfortran do-exit.f90
takk@deb9:~$ ./a.out
           1           1
           2           1
takk@deb9:~$

1列目のIが2回回りましたので、内部のDOだけEXITしたことがわかります。
外部のDOを指定することもできるようです。

takk@deb9:~$ cat do-exit2.f90
      PROGRAM MAIN
        OUTER: DO I=1,2
          DO J=1,3
            IF(J.EQ.2)EXIT OUTER
            WRITE(*,*) I,J
          END DO
        END DO
        STOP
      END
takk@deb9:~$

EXITにOUTERとうラベルを指定して、外部のDOを指定するようにしました。
ビルドしてみます。

takk@deb9:~$ gfortran do-exit2.f90
do-exit2.f90:7:14:

         END DO
              1
Error: Expected block name of ‘outer’ in END DO statement at (1)
do-exit2.f90:9:9:

       END
         1
Error: END DO statement expected at (1)
f951: Error: Unexpected end of file in ‘do-exit2.f90’
takk@deb9:~$

エラーになりました。DOにラベルを付けたら、対応するEND DOにも書かないといけないみたいです。

takk@deb9:~$ cat do-exit2.f90
      PROGRAM MAIN
        OUTER: DO I=1,2
          DO J=1,3
            IF(J.EQ.2)EXIT OUTER
            WRITE(*,*) I,J
          END DO
        END DO OUTER
        STOP
      END
takk@deb9:~$ gfortran do-exit2.f90
takk@deb9:~$

外側のEND DOをEND DO OUTERにしたら、ビルド通りました。
実行してみます。

takk@deb9:~$ ./a.out
           1           1
takk@deb9:~$

外側のDOが指定されたので、一回だけ回って終了しました。
EXITという名前は、他の言語だとプログラム終了だったりするので、紛らわしいかもしれませんが、いろんな言語の兄であるFORTRANが先なんでしょうね。

Fortran90 do cycle

Fortranです。DO CYCLEを使います。
DOの繰り返しを終了するのがEXITでした。前回作ったプログラムのEXITをCYCLEに変えて実行してみます。

takk@deb9:~$ cat do-cycle.f90
      PROGRAM MAIN
        DO I=1,2
          DO J=1,3
            IF(J.EQ.2)CYCLE
            WRITE(*,*) I,J
          END DO
        END DO
        STOP
      END
takk@deb9:~$ gfortran do-cycle.f90
takk@deb9:~$ ./a.out
           1           1
           1           3
           2           1
           2           3
takk@deb9:~$

次は外側のDOにラベルを付けて、CYCLEで指定してみます。

takk@deb9:~$ cat do-cycle2.f90
      PROGRAM MAIN
        OUTER: DO I=1,2
          DO J=1,3
            IF(J.EQ.2)CYCLE OUTER
            WRITE(*,*) I,J
          END DO
        END DO OUTER
        STOP
      END
takk@deb9:~$ gfortran do-cycle2.f90
takk@deb9:~$ ./a.out
           1           1
           2           1
takk@deb9:~$

あれ、この動きはEXITと同じなのでは。

takk@deb9:~$ cat do-exit.f90
      PROGRAM MAIN
        DO I=1,2
          DO J=1,3
            IF(J.EQ.2)EXIT
            WRITE(*,*) I,J
          END DO
        END DO
        STOP
      END
takk@deb9:~$ gfortran do-exit.f90
takk@deb9:~$ ./a.out
           1           1
           2           1
takk@deb9:~$

まったく同じですね。

ところで他の言語はどう書くのでしょう。Pythonで同じようなプログラムを書いてみます。

takk@deb9:~$ cat break.py
for i in range(1,3):
  for j in range(1,4):
    if j == 2:
      break
    print i,j

takk@deb9:~$ python break.py
1 1
2 1
takk@deb9:~$ vi break.py
takk@deb9:~$ cp break.py continue.py
takk@deb9:~$ vi continue.py
takk@deb9:~$ cat continue.py
for i in range(1,3):
  for j in range(1,4):
    if j == 2:
      continue
    print i,j

takk@deb9:~$ python continue.py
1 1
1 3
2 1
2 3
takk@deb9:~$

Fortranのexitが、Pythonのbreakで、Fortranのcycleが、Pythonのcontinueですね。
PythonのこのキーワードはC言語から派生してきたのでしょうか。exitとcycleよりも、breakとcontinueの方が分かりやすい気がします。

Fortran90 大文字小文字区別

Fortranです。Fortran90の学習なので拡張子をf90に変えてビルドしていましたが、プログラム自体を全部大文字で書いてました。小文字で書く方がFortranっぽいらしいので、これからは小文字を使っていこうと思います。ところで90以降は、大文字小文字の区別はあるのでしょうか、実験してみます。

takk@deb9:~$ gfortran var.f90
takk@deb9:~$ cat var.f90
      program main
        I=1
        i=2
        write(*,*) I,i
        stop
      end
takk@deb9:~$ gfortran var.f90
takk@deb9:~$ ./a.out
           2           2
takk@deb9:~$

このように、変数Iとiに別々の値を設定して、表示してみます。

takk@deb9:~$ gfortran var.f90
takk@deb9:~$ ./a.out
           2           2
takk@deb9:~$

同じ値が表示されました。変数名の大文字小文字の区別はないようです。
関数はどうでしょうか。
まずは適当に関数を呼び出しするプログラムを作って、

takk@deb9:~$ cat sub.f90
      program main
        call pr()
        stop
      end
      subroutine pr()
        write(*,*) 'HELLO'
      end
takk@deb9:~$ gfortran sub.f90
takk@deb9:~$ ./a.out
 HELLO
takk@deb9:~$

ビルド、実行確認したら、次にprを大文字にしたPR()を追加します。

takk@deb9:~$ cat sub.f90
      program main
        call pr()
        stop
      end
      subroutine pr()
        write(*,*) 'HELLO'
      end
      subroutine PR()
        write(*,*) 'HELLO'
      end
takk@deb9:~$

そしてビルド。

takk@deb9:~$ gfortran sub.f90
sub.f90:8:19:

sub.f90:5:19:

       subroutine pr()
                   2
sub.f90:8:19:

       subroutine PR()
                   1
Error: Global name ‘pr’ at (1) is already being used as a SUBROUTINE at (2)
takk@deb9:~$

エラーとなりました。やはり大文字小文字の区別はないようです。
gfortranで大文字小文字を区別するオプションは探すことができませんでしたが、他のコンパイラではありそうです。

Fortran90 定数

ずっとFortranは定数が使えないと思ってましたが、ありました。parameterを使うようです。

定数を使う前のプログラム。

takk@deb9:~$ cat do.f90
      program main
        do i=1,2
         do j=1,3
          write(*,*) 'i=',i,' / j=',j
         end do
        end do
        stop
      end
takk@deb9:~$ gfortran do.f90
takk@deb9:~$ ./a.out
 i=           1  / j=           1
 i=           1  / j=           2
 i=           1  / j=           3
 i=           2  / j=           1
 i=           2  / j=           2
 i=           2  / j=           3
takk@deb9:~$

上のプログラムを定数を使った形に変えてみます。

takk@deb9:~$ cat para.f90
      program main
        parameter (ifirst=1)
        parameter (imax=2)
        parameter (jfirst=1)
        parameter (jmax=3)
        do i=ifirst,imax
         do j=jfirst,jmax
          write(*,*) 'i=',i,' / j=',j
         end do
        end do
        stop
      end
takk@deb9:~$ gfortran para.f90
takk@deb9:~$ ./a.out
 i=           1  / j=           1
 i=           1  / j=           2
 i=           1  / j=           3
 i=           2  / j=           1
 i=           2  / j=           2
 i=           2  / j=           3
takk@deb9:~$

変数の代入を()で括って、paramterをつければ定数になるようですが、本当に定数なのでしょうか。代入文を追加してビルドエラーとなるか確認してみます。

takk@deb9:~$ cat para2.f90
      program main
        parameter (ifirst=1)
        parameter (imax=2)
        parameter (jfirst=1)
        parameter (jmax=3)
        imax=5
        do i=ifirst,imax
         do j=jfirst,jmax
          write(*,*) 'i=',i,' / j=',j
         end do
        end do
        stop
      end
takk@deb9:~$ gfortran para2.f90
para2.f90:6:8:

         imax=5
        1
Error: Named constant ‘imax’ in variable definition context (assignment) at (1)
takk@deb9:~$

imax=5という代入を追加しましたが、ビルドでエラーとなりました。ちゃんと定数となっているようですね。

上記までのプログラムでは、parameterを複数行使い複数定義してましたが、変数宣言と同じく、一つのparamterに対して、複数の定数定義もできるのでこのような書き方もできます。

        parameter (ifirst=1,imax=2,jfirst=1,jmax=3)

Fortran90 data

Fortran続きです。dataで配列の初期化を行います。

takk@deb9:~$ cat data.f90
      program main
        integer a(10),s
        data a /5,5,5,5,5,5,5,5,5,5/
        s = 0
        write(*,*) a(3)
        do i = 1,10
          s = s + a(i)
        end do
        write(*,*) s
        stop
      end
takk@deb9:~$
takk@deb9:~$ gfortran data.f90
takk@deb9:~$ ./a.out
           5
          50
takk@deb9:~$

値が同じならまとめて記述することもできます。

takk@deb9:~$ cat data.f90
      program main
        integer a(10),s
        data a /5,5,5,5,5,5,5,5,5,5/
        s = 0
        write(*,*) a(3)
        do i = 1,10
          s = s + a(i)
        end do
        write(*,*) s
        stop
      end
takk@deb9:~$ gfortran data.f90
takk@deb9:~$ ./a.out
           5
          50
takk@deb9:~$

初期値を設定したい配列が複数あるときは、カンマで区切って記述できます。

takk@deb9:~$ cat data.f90
      program main
        integer a(10),b(10),s
        data a,b /10*5,10*2/
        s = 0
        write(*,*) a(3)
        do i = 1,10
          s = s + a(i) + b(i)
        end do
        write(*,*) s
        stop
      end
takk@deb9:~$ gfortran data.f90
takk@deb9:~$ ./a.out
           5
          70
takk@deb9:~$

dataを使わず型宣言で初期化することもできます。

takk@deb9:~$ cat arr.f90
      program main
        integer :: a(10) = (/1,2,3,4,5,6,7,8,9,10/)
        integer s
        s = 0
        write(*,*) a(3)
        do i = 1,10
          s = s + a(i)
        end do
        write(*,*) s
        stop
      end
takk@deb9:~$ gfortran arr.f90
takk@deb9:~$ ./a.out
           3
          55
takk@deb9:~$

Fortran90 type

Fortran続きです。なんとFortranでも構造体が使えるんですね。ということで使ってみます。
typeとend typeの間にメンバ変数を定義します。

takk@deb9:~$ cat type.f90
      program main
        type test
          integer i
          character*10 name
        end type test

        type(test) a
        a = test(1,'aaaa')

        write(*,*) a%name
        write(*,*) a%i
        stop
      end
takk@deb9:~$ gfortran type.f90
takk@deb9:~$ ./a.out
 aaaa
           1
takk@deb9:~$

構造体の変数名%メンバ名でアクセスできます。
当然ですが構造体は型なので、量産できます。

takk@deb9:~$ cat type.f90
      program main
        type test
          integer i
          character*10 name
        end type test

        type(test) a,b,c
        a = test(1,'aaaa')
        b = test(10,'bbbb')
        c = test(100,'cccc')

        write(*,*) a%name
        write(*,*) a%i
        write(*,*) b%name
        write(*,*) b%i
        write(*,*) c%name
        write(*,*) c%i
        stop
      end
takk@deb9:~$ gfortran type.f90
takk@deb9:~$ ./a.out
 aaaa
           1
 bbbb
          10
 cccc
         100
takk@deb9:~$

構造体変数の定義は、途中で追加するのはエラーになるようです。

takk@deb9:~$ cat type2.f90
      program main
        type test
          integer i
          character*10 name
        end type test

        type(test) a
        a = test(1,'aaaa')

        write(*,*) a%name
        write(*,*) a%i

        type(test) b,c
        b = test(10,'bbbb')
        c = test(100,'cccc')

        write(*,*) b%name
        write(*,*) b%i
        write(*,*) c%name
        write(*,*) c%i
        stop
      end
takk@deb9:~$ gfortran type2.f90
type2.f90:13:22:

         type(test) b,c
                      1
Error: Unexpected data declaration statement at (1)
type2.f90:17:21:

         write(*,*) b%name
                     1
Error: Symbol ‘b’ at (1) has no IMPLICIT type
type2.f90:18:21:

         write(*,*) b%i
                     1
Error: Symbol ‘b’ at (1) has no IMPLICIT type
type2.f90:19:21:

         write(*,*) c%name
                     1
Error: Symbol ‘c’ at (1) has no IMPLICIT type
type2.f90:20:21:

         write(*,*) c%i
                     1
Error: Symbol ‘c’ at (1) has no IMPLICIT type
type2.f90:14:12:

         b = test(10,'bbbb')
            1
Error: Can't convert TYPE(test) to REAL(4) at (1)
type2.f90:15:12:

         c = test(100,'cccc')
            1
Error: Can't convert TYPE(test) to REAL(4) at (1)
takk@deb9:~$

まあこれは構造体に限ったことではなく、変数の宣言が最初にないといけないのと同じですね。

takk@deb9:~$ cat var.f90
      program main
        integer i
        i = 1
        write(*,*) i

        integer j
        j = 2
        write(*,*) j
        stop
      end
takk@deb9:~$ gfortran var.f90
var.f90:6:17:

         integer j
                 1
Error: Unexpected data declaration statement at (1)
takk@deb9:~$

Fortran90 where

ずっとFortran学習してますが、なかなか飽きません。飽きるまで続く気がしています。そして新たな発見。where構文。こんな強力な構文があるとは知りませんでした。
今回はwhereを使っていきます。

まずはwhereを使わない例です。doで回して配列aの要素が5より大きかったら、配列bに代入するときに、2倍にするプログラムです。

takk@deb9:~$ cat do.f90
      program main
        integer a(10),b(10)
        data a /1,2,3,4,5,6,7,8,9,10/

        do i=1,10
          if (a(i) > 5) then
            b(i) = a(i) * 2
          else
            b(i) = a(i)
          end if
        end do

        do i=1,10
         write(*,*) b(i)
        end do

        stop
      end
takk@deb9:~$ gfortran do.f90
takk@deb9:~$ ./a.out
           1
           2
           3
           4
           5
          12
          14
          16
          18
          20
takk@deb9:~$

6行目から2倍になってますね。

上のプログラムをwhereを使って作り直すとこうなります。

takk@deb9:~$ cat where.f90
      program main
        integer a(10),b(10)
        data a /1,2,3,4,5,6,7,8,9,10/

        where (a > 5)
         b = a * 2
        else where
         b = a
        end where

        do i=1,10
         write(*,*) b(i)
        end do

        stop
      end
takk@deb9:~$

いやあスゴイ。なんて簡潔に書けるのでしょうか。

実行するともちろん同じ結果になります。

takk@deb9:~$ gfortran where.f90
takk@deb9:~$ ./a.out
           1
           2
           3
           4
           5
          12
          14
          16
          18
          20
takk@deb9:~$

使い方は、whereとend whereの間に配列代入文を挿むだけ。

where 条件
  配列代入文
end where

もしくはこのパターン。

where 条件
  配列代入文1
else where
  配列代入文2
end where

配列代入文以外はだめなのでしょうか。配列代入の代わりにwriteを入れてみます。

takk@deb9:~$ cat where2.f90
      program main
        integer a(10),b(10)
        data a /1,2,3,4,5,6,7,8,9,10/

        where (a > 5)
         write(*,*) 'HELLO'
        end where

        stop
      end
takk@deb9:~$

ビルドエラーとなります。

takk@deb9:~$ gfortran where2.f90
where2.f90:6:27:

          write(*,*) 'HELLO'
                           1
Error: Unexpected WRITE statement in WHERE block at (1)
takk@deb9:~$

もちろんwhere構文内なので配列代入しかできないのであって、if文の中なら問題ないです。

takk@deb9:~$ gfortran if.f90
takk@deb9:~$ ./a.out
 HELLO
takk@deb9:~$ cat if.f90
      program main
        integer a(10)
        data a /1,2,3,4,5,6,7,8,9,10/

        if (a(6) > 5) then
         write(*,*) 'HELLO'
        end if

        stop
      end
takk@deb9:~$ gfortran if.f90
takk@deb9:~$ ./a.out
 HELLO
takk@deb9:~$

普通の代入はどうでしょうか。

takk@deb9:~$ cat where3.f90
      program main
        integer a(10),b(10)
        data a /1,2,3,4,5,6,7,8,9,10/

        where (a > 5)
         i = 1
        end where

        stop
      end
takk@deb9:~$ gfortran where3.f90
where3.f90:6:9:

          i = 1
         1
Error: WHERE assignment target at (1) has inconsistent shape
takk@deb9:~$

whereの中では、普通の代入もエラーとなりました。
やはり配列代入のみしか書くことはできないようです。

Fortran90 where その2

whereを使って配列を初期化する方法を考えてみます。

配列を同じ値で初期化する場合は、たとえば1なら、こんな風に書けると思います。

        where (a .ne. 1)
         a = 1
        end where

whereのあとの a .ne. 1は常に真でもよかったかもしれませんが、書き方が分からないので、1以外なら1で初期化するようにしています。
実行してみます。

takk@deb9:~$ cat where4.f90
      program main
        integer a(10)

        where (a .ne. 1)
         a = 1
        end where

        do i = 1,10
          write(*,*) a(i)
        end do

        stop
      end
takk@deb9:~$ gfortran where4.f90
takk@deb9:~$ ./a.out
           1
           1
           1
           1
           1
           1
           1
           1
           1
           1
takk@deb9:~$

1で初期化できました。

data文で初期化するのとどちらがプログラムサイズが少なくて済むでしょうか。配列数を10000ぐらいにして確認してみます。
まずはdata文を使って静的に初期化した場合。

takk@deb9:~$ cat data5.f90
      program main
        integer a(10000)
        data a /10000*1/

        write(*,*) a(5000)
        stop
      end
takk@deb9:~$ gfortran data5.f90
takk@deb9:~$ ./a.out
           1
takk@deb9:~$ ls -l a.out
-rwxr-xr-x 1 takk takk 49136 12月 20 20:34 a.out
takk@deb9:~$

簡潔に書けますね。サイズは49136。大きいです。

where文で動的に初期化した場合。

takk@deb9:~$ cat where5.f90
      program main
        integer a(10000)

        where (a .ne. 1)
         a = 1
        end where

        write(*,*) a(5000)

        stop
      end
takk@deb9:~$ gfortran where5.f90
takk@deb9:~$ ./a.out
           1
takk@deb9:~$ ls -l a.out
-rwxr-xr-x 1 takk takk 9104 12月 20 20:33 a.out
takk@deb9:~$

行数は2行増えましたが、こちらも簡潔に書けてる気がします。プログラムサイズは9104。5分の1になりました。

Fortran90 where その3

前回whereを使って、配列に初期値を格納したのですが、whereを使う意味はまったくありませんでした。
配列代入の意味を勘違いしていたのですが、まずは前回作った初期化のプログラムを見てみます。

takk@deb9:~$ cat where4.f90
      program main
        integer a(10)

        where (a .ne. 1)
         a = 1
        end where

        do i = 1,10
          write(*,*) a(i)
        end do

        stop
      end
takk@deb9:~$

10個の要素すべてが1に初期化されていたので正しく動いてはいたのですが、実はwhereを使う必要はありませんでした。こうすれば良いです。

takk@deb9:~$ cat where6.f90
      program main
        integer a(10)

        a = 1

        do i=1,10
          write(*,*) a(i)
        end do

        stop
      end
takk@deb9:~$

そうです。a = 1とするだけ。10個の要素すべてが1になります。
実行してみます。

takk@deb9:~$ gfortran where6.f90
takk@deb9:~$ ./a.out
           1
           1
           1
           1
           1
           1
           1
           1
           1
           1
takk@deb9:~$

whereはすごい! と糠喜びしてましたが、こうなると、whereって何に使うのだろう? と疑問符。

要素の総和とかできるでしょうか。1~10までの足し算を試みます。

takk@deb9:~$ cat where7.f90
      program main
        integer a(10),s
        data a /1,2,3,4,5,6,7,8,9,10/

        s = 0

        where(a > 0)
          s = s + a
        end where

        do i=1,10
          write(*,*) a(i)
        end do

        stop
      end
takk@deb9:~$

こんなんでいけるでしょうか。ビルドしてみます。

takk@deb9:~$ gfortran where7.f90
where7.f90:8:10:

           s = s + a
          1
Error: Incompatible ranks 0 and 1 in assignment at (1)
where7.f90:8:10:

           s = s + a
          1
Error: WHERE assignment target at (1) has inconsistent shape
takk@deb9:~$

エラーになりました。スカラー変数への代入はwhere構文内ではできないことを忘れてました。
配列代入しかできないということは、総和もできないってことですね。

Fortran90 多次元配列

Fortranで多次元配列を使ってみます。まずは二次元から。

takk@deb9:~$ cat arr1.f90
      program main
        integer a(2,3)
        data a /10,20,30,40,50,60/

        write(*,*) a(1,1)
        write(*,*) a(1,2)
        write(*,*) a(1,3)
        write(*,*) a(2,2)
        write(*,*) a(2,1)
        write(*,*) a(2,3)

        stop
      end
takk@deb9:~$

dataを使って、初期化をしていますが、どのような順番で格納されたのでしょうか。

takk@deb9:~$ gfortran arr1.f90
takk@deb9:~$ ./a.out
          10
          30
          50
          40
          20
          60
takk@deb9:~$

実行結果から、このように格納されているようです。

 (,1) (,2) (,3)
+----+----+----+
| 10 | 30 | 50 | (1,)
+----+----+----+
| 20 | 40 | 60 | (2,)
+----+----+----+

二次元配列も配列代入が使えるのでしょうか。試してみます。

takk@deb9:~$ cat arr2.f90
      program main
        integer a(2,3)

        a = 33

        do m=1,3
          do n=1,2
            write(*,*) a(n,m)
          end do
        end do

        stop
      end
takk@deb9:~$ gfortran arr2.f90
takk@deb9:~$ ./a.out
          33
          33
          33
          33
          33
          33
takk@deb9:~$

できました。1行の代入文ですべて初期化されました。

Fortran90 多次元配列 その2

配列代入が中で何をしているか気になって仕方がありません。GDBを使って見ていきます。
まずは、配列要素が1つで、次元数1つのプログラムで確認してみます。

a(1)

という整数を定義し、配列代入a=100を実行するプログラムです。

takk@deb9:~$ cat arr1x1.f90
      program main
        integer a(1)
        a = 100
        stop
      end
takk@deb9:~$ gfortran -g arr1x1.f90
takk@deb9:~$

GDBを起動します。

takk@deb9:~$ gfortran -g arr1x1.f90
takk@deb9:~$ gdb a.out

          ~

(gdb) l
1             program main
2               integer a(1)
3               a=100
4               stop
5             end
(gdb)

a=100の行にブレイクポイントを設定し実行します。

(gdb) b 3
Breakpoint 1 at 0x7d8: file arr1x1.f90, line 3.
(gdb) r
Starting program: /home/takk/a.out

Breakpoint 1, MAIN__ () at arr1x1.f90:3
3               a=100
(gdb) 

アセンブラを見てみましょう。

(gdb) disas
Dump of assembler code for function MAIN__:
   0x00005555555547d0 <+0>:     push   %rbp
   0x00005555555547d1 <+1>:     mov    %rsp,%rbp
   0x00005555555547d4 <+4>:     sub    $0x10,%rsp
=> 0x00005555555547d8 <+8>:     mov    $0x1,%eax
   0x00005555555547dd <+13>:    cmp    $0x1,%rax
   0x00005555555547e1 <+17>:    jg     0x5555555547f5 <MAIN__+37>
   0x00005555555547e3 <+19>:    lea    -0x1(%rax),%rdx
   0x00005555555547e7 <+23>:    movl   $0x64,-0x4(%rbp,%rdx,4)
   0x00005555555547ef <+31>:    add    $0x1,%rax
   0x00005555555547f3 <+35>:    jmp    0x5555555547dd <MAIN__+13>
   0x00005555555547f5 <+37>:    mov    $0x0,%esi
   0x00005555555547fa <+42>:    mov    $0x0,%edi
   0x00005555555547ff <+47>:    callq  0x555555554660 <_gfortran_stop_string@plt>
End of assembler dump.
(gdb)

100は、16進で0x64ですので、<+23>の行あたりで、a(1)へ配列代入しているものと思われます。<+35>の行のjmpは、アドレス0x5555555547ddに戻っているので配列のループを作っていますね。同アドレスでチェックしてます。
この作りを覚えておいて、次は二次元配列にした場合、どう変わるのか見てみます。

takk@deb9:~$ cat arr1x2.f90
      program main
        integer a(1,1)
        a=100
        stop
      end
takk@deb9:~$

アセンブラを見てみます。

(gdb) disas
Dump of assembler code for function MAIN__:
   0x00005555555547d0 <+0>:     push   %rbp
   0x00005555555547d1 <+1>:     mov    %rsp,%rbp
   0x00005555555547d4 <+4>:     sub    $0x10,%rsp
=> 0x00005555555547d8 <+8>:     mov    $0x1,%eax
   0x00005555555547dd <+13>:    cmp    $0x1,%rax
   0x00005555555547e1 <+17>:    jg     0x55555555480a <MAIN__+58>
   0x00005555555547e3 <+19>:    lea    -0x2(%rax),%rsi
   0x00005555555547e7 <+23>:    mov    $0x1,%edx
   0x00005555555547ec <+28>:    cmp    $0x1,%rdx
   0x00005555555547f0 <+32>:    jg     0x555555554804 <MAIN__+52>
   0x00005555555547f2 <+34>:    lea    (%rdx,%rsi,1),%rcx
   0x00005555555547f6 <+38>:    movl   $0x64,-0x4(%rbp,%rcx,4)
   0x00005555555547fe <+46>:    add    $0x1,%rdx
   0x0000555555554802 <+50>:    jmp    0x5555555547ec <MAIN__+28>
   0x0000555555554804 <+52>:    add    $0x1,%rax
   0x0000555555554808 <+56>:    jmp    0x5555555547dd <MAIN__+13>
   0x000055555555480a <+58>:    mov    $0x0,%esi
   0x000055555555480f <+63>:    mov    $0x0,%edi
   0x0000555555554814 <+68>:    callq  0x555555554660 <_gfortran_stop_string@plt>
End of assembler dump.
(gdb)

前に戻るjmpが2つになったので、次元数と一致します。要素数が1でもそれ以外でも関係なくこの作りになってそうです。
ではビルドの限界である七次元配列で確認してみます。要素数はそれぞれ2にしておきます。
確認に使ったプログラムはこれです。

takk@deb9:~$ cat arr2x7.f90
      program main
        integer a(2,2,2,2,2,2,2)
        a=100
        stop
      end
takk@deb9:~$

アセンブラです。

 (gdb) disas
Dump of assembler code for function MAIN__:
   0x00005555555547d0 <+0>:     push   %rbp
   0x00005555555547d1 <+1>:     mov    %rsp,%rbp
   0x00005555555547d4 <+4>:     push   %r15
   0x00005555555547d6 <+6>:     push   %r14
   0x00005555555547d8 <+8>:     push   %r13
   0x00005555555547da <+10>:    push   %r12
   0x00005555555547dc <+12>:    push   %rbx
   0x00005555555547dd <+13>:    sub    $0x208,%rsp
=> 0x00005555555547e4 <+20>:    mov    $0x1,%edi
   0x00005555555547e9 <+25>:    cmp    $0x2,%rdi
   0x00005555555547ed <+29>:    jg     0x5555555548c3 <MAIN__+243>
   0x00005555555547f3 <+35>:    mov    %rdi,%rax
   0x00005555555547f6 <+38>:    shl    $0x6,%rax
   0x00005555555547fa <+42>:    lea    -0x7f(%rax),%rbx
   0x00005555555547fe <+46>:    mov    $0x1,%r8d
   0x0000555555554804 <+52>:    cmp    $0x2,%r8
   0x0000555555554808 <+56>:    jg     0x5555555548ba <MAIN__+234>
   0x000055555555480e <+62>:    mov    %r8,%rax
   0x0000555555554811 <+65>:    shl    $0x5,%rax
   0x0000555555554815 <+69>:    lea    (%rax,%rbx,1),%r12
   0x0000555555554819 <+73>:    mov    $0x1,%r9d
   0x000055555555481f <+79>:    cmp    $0x2,%r9
   0x0000555555554823 <+83>:    jg     0x5555555548b1 <MAIN__+225>
   0x0000555555554829 <+89>:    mov    %r9,%rax
   0x000055555555482c <+92>:    shl    $0x4,%rax
   0x0000555555554830 <+96>:    lea    (%rax,%r12,1),%r13
   0x0000555555554834 <+100>:   mov    $0x1,%eax
   0x0000555555554839 <+105>:   cmp    $0x2,%rax
   0x000055555555483d <+109>:   jg     0x5555555548a8 <MAIN__+216>
   0x000055555555483f <+111>:   lea    0x0(,%rax,8),%rdx
   0x0000555555554847 <+119>:   lea    (%rdx,%r13,1),%r14
   0x000055555555484b <+123>:   mov    $0x1,%edx
   0x0000555555554850 <+128>:   cmp    $0x2,%rdx
   0x0000555555554854 <+132>:   jg     0x5555555548a2 <MAIN__+210>
   0x0000555555554856 <+134>:   lea    0x0(,%rdx,4),%rcx
   0x000055555555485e <+142>:   lea    (%rcx,%r14,1),%r15
   0x0000555555554862 <+146>:   mov    $0x1,%ecx
   0x0000555555554867 <+151>:   cmp    $0x2,%rcx
   0x000055555555486b <+155>:   jg     0x55555555489c <MAIN__+204>
   0x000055555555486d <+157>:   lea    (%rcx,%rcx,1),%rsi
   0x0000555555554871 <+161>:   lea    (%rsi,%r15,1),%r11
   0x0000555555554875 <+165>:   mov    $0x1,%esi
   0x000055555555487a <+170>:   cmp    $0x2,%rsi
   0x000055555555487e <+174>:   jg     0x555555554896 <MAIN__+198>
   0x0000555555554880 <+176>:   lea    (%rsi,%r11,1),%r10
   0x0000555555554884 <+180>:   movl   $0x64,-0x230(%rbp,%r10,4)
   0x0000555555554890 <+192>:   add    $0x1,%rsi
   0x0000555555554894 <+196>:   jmp    0x55555555487a <MAIN__+170>
   0x0000555555554896 <+198>:   add    $0x1,%rcx
   0x000055555555489a <+202>:   jmp    0x555555554867 <MAIN__+151>
   0x000055555555489c <+204>:   add    $0x1,%rdx
   0x00005555555548a0 <+208>:   jmp    0x555555554850 <MAIN__+128>
   0x00005555555548a2 <+210>:   add    $0x1,%rax
   0x00005555555548a6 <+214>:   jmp    0x555555554839 <MAIN__+105>
   0x00005555555548a8 <+216>:   add    $0x1,%r9
   0x00005555555548ac <+220>:   jmpq   0x55555555481f <MAIN__+79>
   0x00005555555548b1 <+225>:   add    $0x1,%r8
   0x00005555555548b5 <+229>:   jmpq   0x555555554804 <MAIN__+52>
   0x00005555555548ba <+234>:   add    $0x1,%rdi
   0x00005555555548be <+238>:   jmpq   0x5555555547e9 <MAIN__+25>
   0x00005555555548c3 <+243>:   mov    $0x0,%esi
   0x00005555555548c8 <+248>:   mov    $0x0,%edi
   0x00005555555548cd <+253>:   callq  0x555555554660 <_gfortran_stop_string@plt>
End of assembler dump.
(gdb)

予想通り配列処理のループを作っているjmpが7つあります。(jmpqが混ざってますが)
配列代入といえども、一つずつ代入していることには変わりありませんね。

Fortran90 多次元配列 その3

配列の値を表示するのにwriteを使うため、doで繰り返し制御していましたが、簡単な書き方があるようです。使ってみます。
まずは、今までのdoでwriteを繰り返し実行する方法。

takk@deb9:~$ cat arr1.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        do i=1,5
          write(*,*) a(i)
        end do

        stop
      end
takk@deb9:~$ gfortran arr1.f90
takk@deb9:~$ ./a.out
           1
           2
           3
           4
           5
takk@deb9:~$

次にdoを使わず、writeに指定する変数に指定する方法。do型反復というらしいです。

a(i)の部分を、
(a(i),i=1,5)に置き換えます。

takk@deb9:~$ cat arr2.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,*) (a(i),i=1,5)

        stop
      end
takk@deb9:~$ gfortran arr2.f90
takk@deb9:~$ ./a.out
           1           2           3           4           5
takk@deb9:~$

タブ区切りで全要素表示されました。
これはこれで良いのですが、改行して表示したい場合どうすれば良いでしょう。

takk@deb9:~$ cat arr3.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,10010) (a(i),i=1,5)

        stop
10010 format(I3)
      end
takk@deb9:~$ gfortran arr3.f90
takk@deb9:~$ ./a.out
  1
  2
  3
  4
  5
takk@deb9:~$

formatで指定するだけでした。

Fortran90 行列

行列計算します。まずは行列の和から。
配列代入を使えば簡単です。

takk@deb9:~$ cat matrix-sum.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /1,2, 3,4/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)

        c = a + b

        write(*,*) 'c=a+b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$

引き算も同じ要領です。

takk@deb9:~$ cat matrix-sub.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /11,22, 33,44/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)
        c = a - b

        write(*,*) 'c=a-b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$ gfortran matrix-sub.f90
takk@deb9:~$ ./a.out
 a=
          11          33
          22          44
 b=
          10          30
          20          40
 c=a-b
           1           3
           2           4
takk@deb9:~$

では積はどうでしょう。a * bだけで済むのでしょうか。

takk@deb9:~$ cat matrix-mul.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /1,2, 3,4/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)
        c = a * b

        write(*,*) 'c=a*b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$ gfortran matrix-mul.f90
takk@deb9:~$ ./a.out
 a=
           1           3
           2           4
 b=
          10          30
          20          40
 c=a*b
          10          90
          40         160
takk@deb9:~$

配列の要素どうしの積になってしまいました。
このような計算をしてほしかったのですが、配列代入だけではできないようです。
\[
\left(
\begin{array}{cc}
a1 & a2 \\
a3 & a4
\end{array}
\right)
\times \left(
\begin{array}{cc}
b1 & b2 \\
b3 & b4
\end{array}
\right)
= \left(
\begin{array}{cc}
a1 \times b1 + a2 \times b3 & a1 \times b2 + a2 \times b4 \\
a3 \times b1 + a4 \times b3 & a3 \times b2 + a4 \times b4
\end{array}
\right)
\]

次回は、この行列の積を作ってみます。

Fortran90 行列の積

前回は行列の内積になってしまったので、通常の積をしてみます。この式でしたね。
\[
\left(
\begin{array}{cc}
a11 & a12 \\
a21 & a22
\end{array}
\right)
\times \left(
\begin{array}{cc}
b11 & b12 \\
b21 & b22
\end{array}
\right)
= \left(
\begin{array}{cc}
a11 \times b11 + a12 \times b21 & a11 \times b12 + a12 \times b22 \\
a21 \times b11 + a22 \times b21 & a21 \times b12 + a22 \times b22
\end{array}
\right)
\]
式をそのままプログラムしてみます。

takk@deb9:~$ cat matrix-mul.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /1,2, 3,4/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)

        c(1,1) = a(1,1) * b(1,1) + a(1,2) * b(2,1)
        c(1,2) = a(1,1) * b(1,2) + a(1,2) * b(2,2)
        c(2,1) = a(2,1) * b(1,1) + a(2,2) * b(2,1)
        c(2,2) = a(2,1) * b(1,2) + a(2,2) * b(2,2)

        write(*,*) 'c=a*b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$

実行すると、

takk@deb9:~$ gfortran matrix-mul.f90
takk@deb9:~$ ./a.out
 a=
           1           3
           2           4
 b=
          10          30
          20          40
 c=a*b
          70         150
         100         220
takk@deb9:~$

行列の積をするmatmulという組み込み関数があるらしいので、使ってみます。

takk@deb9:~$ cat matrix-mul2.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /1,2, 3,4/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)

        c = matmul(a,b)

        write(*,*) 'c=a*b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$

実行すると、

takk@deb9:~$ gfortran matrix-mul2.f90
takk@deb9:~$ ./a.out
 a=
           1           3
           2           4
 b=
          10          30
          20          40
 c=a*b
          70         150
         100         220
takk@deb9:~$

c=a*bの答えが同じになりました。

Fortran90 sum

前回使ったmatmulは配列で使える関数でしたが、他にも配列用関数はあります。
今回は、sumを使ってみます。
まずはsumを使わないプログラムです。

takk@deb9:~$ cat func-sum1.f90
      program main
        integer a(5),s
        data a /1,2,3,4,5/

        s = 0
        do i = 1,5
          s = s + a(i)
        end do
        write(*,*) s

        stop
      end
takk@deb9:~$ gfortran func-sum1.f90
takk@deb9:~$ ./a.out
          15
takk@deb9:~$

1~5までの和なので結果は15で合ってますね。

takk@deb9:~$ cat func-sum2.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,*) sum(a)

        stop
      end
takk@deb9:~$

ものすごく簡単になりました。
結果はこれです。

takk@deb9:~$ gfortran func-sum2.f90
takk@deb9:~$ ./a.out
          15
takk@deb9:~$

配列の要素の総和ですので、初期化されていない要素があるとこういうことになります。

takk@deb9:~$ cat func-sum3.f90
      program main
        integer a(10)

        a(1)=1
        a(2)=2
        a(3)=3
        a(4)=4
        a(5)=5

        write(*,*) sum(a)

        stop
      end
takk@deb9:~$ gfortran func-sum3.f90
takk@deb9:~$ ./a.out
       65442
takk@deb9:~$

初期化されていない値と計算されて、よく分からない結果になりました。
このようなときは、計算したい要素の範囲を指定すればよいです。

takk@deb9:~$ cat func-sum4.f90
      program main
        integer a(10)

        a(1)=1
        a(2)=2
        a(3)=3
        a(4)=4
        a(5)=5

        write(*,*) sum(a(1:5))

        stop
      end
takk@deb9:~$ gfortran func-sum4.f90
takk@deb9:~$ ./a.out
          15
takk@deb9:~$

Fortran90 sum その2

多次元配列のsumをしてみます。このように値を格納した2×5の配列の要素の総和を計算します。

+--+--+--+--+--+
| 1| 2| 3| 4| 5|
+--+--+--+--+--+
|10|20|30|40|50|
+--+--+--+--+--+

sumの引数に配列名を指定するだけで求めることができます。

takk@deb9:~$ cat func-sum-2-1.f90
      program main
        integer a(2,5)
        data a /1,10,2,20,3,30,4,40,5,50/

        write(*,*) sum(a)

        stop
      end
takk@deb9:~$
takk@deb9:~$ gfortran func-sum-2-1.f90
takk@deb9:~$ ./a.out
         165
takk@deb9:~$

dim=で指定することで、次元別の計算もできます。

takk@deb9:~$ cat func-sum-2-2.f90
      program main
        integer a(2,5)
        data a /1,10,2,20,3,30,4,40,5,50/

        write(*,*) sum(a,dim=1)
        write(*,*) sum(a,dim=2)

        stop
      end
takk@deb9:~$
takk@deb9:~$ gfortran func-sum-2-2.f90
takk@deb9:~$ ./a.out
          11          22          33          44          55
          15         150
takk@deb9:~$

mask=を使うことで、各要素のフィルタリングもすることができます。
10未満の値だけsumするプログラムです。

takk@deb9:~$ cat func-sum-2-3.f90
      program main
        integer a(2,5)
        data a /1,10,2,20,3,30,4,40,5,50/

        write(*,*)sum(a,mask=a .lt. 10)

        stop
      end
takk@deb9:~$

maskを使って偶数の要素の総和、奇数の要素の総和をそれぞれ求めてみます。

takk@deb9:~$ cat func-sum-2-4.f90
      program main
        integer a(10)
        data a /1,2,3,4,5,6,7,8,9,10/

        write(*,*)sum(a,mask=mod(a,2) .eq. 0)
        write(*,*)sum(a,mask=mod(a,2) .eq. 1)

        stop
      end
takk@deb9:~$

結果です。

takk@deb9:~$ gfortran func-sum-2-4.f90
takk@deb9:~$ ./a.out
          30
          25
takk@deb9:~$

Fortran90 product

productを使って配列の積を求めます。
全要素を掛けた結果、つまり総乗を求めたいので、このような計算となります。

takk@deb9:~$ cat func-product-1.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,*) a(1) * a(2) * a(3) * a(4) * a(5)

        stop
      end
takk@deb9:~$ gfortran func-product-1.f90
takk@deb9:~$ ./a.out
         120
takk@deb9:~$

productを使えば一発です。

takk@deb9:~$ cat func-product-2.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,*) product(a)

        stop
      end
takk@deb9:~$ gfortran func-product-2.f90
takk@deb9:~$ ./a.out
         120
takk@deb9:~$

sumと同じくdimやmaskを指定してフィルタリングすることができます。

takk@deb9:~$ cat func-product-3.f90
      program main
        integer a(2,5)
        data a /1,10,2,20,3,30,4,40,5,50/

        write(*,*) product(a,dim=1)
        write(*,*) product(a,dim=2,mask=a.le.30)

        stop
      end
takk@deb9:~$ gfortran func-product-3.f90
takk@deb9:~$ ./a.out
          10          40          90         160         250
         120        6000
takk@deb9:~$

product関数というと、他言語だと掛け合わせのリストを生成するタイプだと思いますが、Fortranでは単純に掛け算なんですね。sum関数の総和に対して、product関数の総乗なので、こちらの関数仕様の方が分かりやすい気がします。

Fortran90 shape,maxval,minval,reshape

配列で使えるその他の関数としてshape,maxval,minval,reshapeを使ってみます。
まずはshapeから。shapeは、配列の形状を戻します。

takk@deb9:~$ cat func-shape.f90
      program main
        integer a(5)
        integer b(2,3)

        write(*,*) shape(a)
        write(*,*) shape(b)

        b = 0
        write(*,*) shape(b)

        stop
      end
takk@deb9:~$

shapeは値が未初期化でも初期化していても、常に配列の各次元の要素数を戻してくれます。

takk@deb9:~$ gfortran func-shape.f90
takk@deb9:~$ ./a.out
           5
           2           3
           2           3
takk@deb9:~$

次は要素のMAX値、MIN値を調べる関数を使います。sumやproductと同じく、dim指定もmask指定もできて便利です。

takk@deb9:~$ cat func-val.f90
      program main
        integer a(2,5)
        data a /1,10,2,20,3,30,4,40,5,50/

        write(*,*) 'max'
        write(*,*) maxval(a)
        write(*,*) maxval(a,dim=1)
        write(*,*) maxval(a,dim=2)

        write(*,*) 'min'
        write(*,*) minval(a)
        write(*,*) minval(a,dim=1)
        write(*,*) minval(a,dim=2)

        write(*,*) 'mask'
        write(*,*) maxval(a,mask=a.ne.50)
        write(*,*) minval(a,mask=a.ne.1)

        stop
      end
takk@deb9:~$

結果です。

takk@deb9:~$ gfortran func-val.f90
takk@deb9:~$ ./a.out
 max
          50
          10          20          30          40          50
           5          50
 min
           1
           1           2           3           4           5
           1          10
 mask
          40
           2
takk@deb9:~$

次は配列を別の形状の配列に置き換えるreshapeです。使いこなせるようになるとかなり自由にプログラミングできるようになるのではないか、と思います。

takk@deb9:~$ cat func-reshape.f90
      program main
        integer a(6),b(2,3)
        data a /1,2,3,4,5,6/

        b = reshape(a,(/2,3/))

        write(*,*) b(1,3)

        stop
      end
takk@deb9:~$
takk@deb9:~$ gfortran func-reshape.f90
takk@deb9:~$ ./a.out
           5
takk@deb9:~$

Fortran90 cshift,eoshift

配列の要素のシフトをしてくれる関数cshiftとeoshiftを使います。
cshiftはcircular shiftを意味しており、循環するシフトです。

takk@deb9:~$ cat func-cshift.f90
      program main
        integer a(5),b(5)
        data a /1,2,3,4,5/

        write(*,*) a
        b = cshift(a,2)
        write(*,*) b

        stop
      end
takk@deb9:~$ gfortran func-cshift.f90
takk@deb9:~$ ./a.out
           1           2           3           4           5
           3           4           5           1           2
takk@deb9:~$

上記ではshift値に2を指定したので要素が2つ左にずれて、飛び出した2個の要素は、右側に収まりました。

shift値をマイナスにすると逆方向にシフトします。

takk@deb9:~$ cat func-cshift2.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,*) 'left'
        write(*,*) cshift(a,1)
        write(*,*) cshift(a,2)
        write(*,*) cshift(a,3)
        write(*,*) cshift(a,4)
        write(*,*) 'right'
        write(*,*) cshift(a,-1)
        write(*,*) cshift(a,-2)
        write(*,*) cshift(a,-3)
        write(*,*) cshift(a,-4)

        stop
      end
takk@deb9:~$
takk@deb9:~$ gfortran func-cshift2.f90
takk@deb9:~$ ./a.out
 left
           2           3           4           5           1
           3           4           5           1           2
           4           5           1           2           3
           5           1           2           3           4
 right
           5           1           2           3           4
           4           5           1           2           3
           3           4           5           1           2
           2           3           4           5           1
takk@deb9:~$

eoshiftは、end-off shiftの略だそうです。こちらは循環しません。シフトで空になった場所には、値を指定しないと0が入ります。

takk@deb9:~$ cat func-eoshift.f90
      program main
        integer a(5),b(5)
        data a /1,2,3,4,5/

        write(*,*) a
        b = eoshift(a,2)
        write(*,*) b

        stop
      end
takk@deb9:~$ gfortran func-eoshift.f90
takk@deb9:~$ ./a.out
           1           2           3           4           5
           3           4           5           0           0
takk@deb9:~$

シフトで空になった要素に値を指定するには、boundary=を指定します。

takk@deb9:~$ cat func-eoshift2.f90
      program main
        integer a(5)
        data a /1,2,3,4,5/

        write(*,*) eoshift(a,1,boundary=999)
        write(*,*) eoshift(a,2,boundary=999)
        write(*,*) eoshift(a,3,boundary=999)
        write(*,*) eoshift(a,4,boundary=999)
        write(*,*) eoshift(a,-1,boundary=999)
        write(*,*) eoshift(a,-2,boundary=999)
        write(*,*) eoshift(a,-3,boundary=999)
        write(*,*) eoshift(a,-4,boundary=999)

        stop
      end
takk@deb9:~$ gfortran func-eoshift2.f90
takk@deb9:~$ ./a.out
           2           3           4           5         999
           3           4           5         999         999
           4           5         999         999         999
           5         999         999         999         999
         999           1           2           3           4
         999         999           1           2           3
         999         999         999           1           2
         999         999         999         999           1
takk@deb9:~$

Fortran90 merge

mergeという言葉をずっと誤解して覚えていたので、このmerge関数も単に配列を結合するだけと思っていました。
おそらくこんな使い方だと思ったので、

aとbの結合後の配列 = merge(a,b)

下のようにプログラミングしてビルドしようとすると、

takk@deb9:~$ cat func-merge-1.f90
      program main
        integer a(3),b(3)
        data a /1,2,3/
        data b /4,5,6/

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

        stop
      end
takk@deb9:~$

ビルドエラーになりました。関数の引数(mask)がないようです。

takk@deb9:~$ gfortran func-merge-1.f90
func-merge-1.f90:8:28:

         write(*,*) 'merge=',merge(a,b)
                            1
Error: Missing actual argument ‘mask’ in call to ‘merge’ at (1)
takk@deb9:~$

とりあえず全部.true.を指定してビルド、実行してみましょう。

takk@deb9:~$ cat func-merge-2.f90
      program main
        integer a(3),b(3)
        data a /1,2,3/
        data b /4,5,6/

        write(*,*) 'a=',a
        write(*,*) 'b=',b
        write(*,*) 'merge=',merge(a,b,mask=.true.)

        stop
      end
takk@deb9:~$
takk@deb9:~$ gfortran func-merge-2.f90
takk@deb9:~$ ./a.out
 a=           1           2           3
 b=           4           5           6
 merge=           1           2           3
takk@deb9:~$

おかしいです。merge=の行に表示されている配列がaの配列そのままです。
それもそのはず。マージは結合(join)ではなく、併合なので、配列サイズが変わるわけではなかったんです。どちらの配列の値を採用するかmaskで決めて配列を作り直す関数なんですね。
このように使います。

takk@deb9:~$ cat func-merge-3.f90
      program main
        integer a(3),b(3)
        data a /10,2,30/
        data b /4,50,6/

        write(*,*) 'a=',a
        write(*,*) 'b=',b
        write(*,*) 'merge=',merge(a,b,mask=a>b)

        stop
      end
takk@deb9:~$

mask=a>bなので、aの要素の方がbより大きければ、aの要素を使い、そうでなければ、bの要素を使う指定となります。要するに大きい方の値でマージするということです。
実行してみます。

takk@deb9:~$ gfortran func-merge-3.f90
takk@deb9:~$ ./a.out
 a=          10           2          30
 b=           4          50           6
 merge=          10          50          30
takk@deb9:~$

Fortran90 spread

spread関数を使います。
spreadは、配列の次元を増やす関数です。増えた配列の値は、元の配列からコピーされます。

spread(配列, 次元,コピー数)
takk@deb9:~$ cat func-spread-1.f90
      program main
        integer a(3)
        data a /1,2,3/

        write(*,*) a
        write(*,*) spread(a,dim=2,ncopies=5)

        stop
      end
takk@deb9:~$ gfortran func-spread-1.f90
takk@deb9:~$ ./a.out
           1           2           3
           1           2           3           1           2           3           1           2           3           1           2           3           1           2           3
takk@deb9:~$

結果がよくわからないので、shape関数を使ってspreadした後の配列の形状を確認します。

takk@deb9:~$ cat func-spread-2.f90
      program main
        integer a(3)
        data a /1,2,3/

        write(*,*) shape(spread(a,dim=2,ncopies=5))

        stop
      end
takk@deb9:~$ gfortran func-spread-2.f90
takk@deb9:~$ ./a.out
           3           5
takk@deb9:~$

3行5列の配列となりました。このようにコピーされた値が格納された配列となりました。

+-+-+-+-+-+
|1|1|1|1|1|
+-+-+-+-+-+
|2|2|2|2|2|
+-+-+-+-+-+
|3|3|3|3|3|
+-+-+-+-+-+

dim=を1に変えると、配列を拡張する方向(次元)が1に変わります。

takk@deb9:~$ cat func-spread-3.f90
      program main
        integer a(3)
        data a /1,2,3/

        write(*,*) shape(spread(a,dim=1,ncopies=5))

        stop
      end
takk@deb9:~$ gfortran func-spread-3.f90
takk@deb9:~$ ./a.out
           5           3
takk@deb9:~$

配列のイメージです。

+-+-+-+
|1|2|3|
+-+-+-+
|1|2|3|
+-+-+-+
|1|2|3|
+-+-+-+
|1|2|3|
+-+-+-+
|1|2|3|
+-+-+-+

do文を使って、イメージ通り格納されているか確認しみましょう。

takk@deb9:~$ cat func-spread-4.f90
      program main
        integer a(3),b(3,5),c(5,3)
        data a /1,2,3/

        b = spread(a,dim=2,ncopies=5)
        c = spread(a,dim=1,ncopies=5)

        write(*,*) 'b(3,5)'
        do i=1,3
            write(*,*) (b(i,j),j=1,5)
        end do

        write(*,*) 'c(5,3)'
        do i=1,5
            write(*,*) (c(i,j),j=1,3)
        end do

        stop
      end
takk@deb9:~$ gfortran func-spread-4.f90
takk@deb9:~$ ./a.out
 b(3,5)
           1           1           1           1           1
           2           2           2           2           2
           3           3           3           3           3
 c(5,3)
           1           2           3
           1           2           3
           1           2           3
           1           2           3
           1           2           3
takk@deb9:~$

Fortran90 transpose

Fortranでtransposeを使います。transposeは、Rubyにもありますね。Fortranを参考にしたのでしょうか。
Ruby(irb)でtransposeを使ってみます。

takk@deb9:~$ irb
irb(main):001:0> a=Array.new(2,Array.new(3))
=> [[nil, nil, nil], [nil, nil, nil]]
irb(main):002:0> a.transpose()
=> [[nil, nil], [nil, nil], [nil, nil]]
irb(main):003:0>

2×3の配列が、3×2の配列に転換されました。引数がいらないので簡単に使えますね。
次はFortranのtransposeを使ってみます。

takk@deb9:~$ cat func-transpose.f90
      program main
        integer a(2,5)

        write(*,*) shape(a)
        write(*,*) shape(transpose(a))

        stop
      end
takk@deb9:~$

配列の形状を比較して確認したいのでshapeを表示するプログラムにしました。
実行してみます。

takk@deb9:~$ gfortran func-transpose.f90
takk@deb9:~$ ./a.out
           2           5
           5           2
takk@deb9:~$

2×5の配列が、5×2になりました。Rubyと同じく簡単でしたね。
次に気になるのが、次元数を増やすとどうなるか。三次元の配列で確認してみましょう。

takk@deb9:~$ cat func-transpose-2.f90
      program main
        integer a(2,3,4)

        write(*,*) shape(a)
        write(*,*) shape(transpose(a))

        stop
      end
takk@deb9:~$

ビルドします。

takk@deb9:~$ gfortran func-transpose-2.f90
func-transpose-2.f90:5:35:

         write(*,*) shape(transpose(a))
                                   1
Error: ‘matrix’ argument of ‘transpose’ intrinsic at (1) must be of rank 2
takk@deb9:~$

ビルドエラーになりました。二次元配列が対象なんですね。
Rubyの方はどうでしょう。

irb(main):001:0> a=Array.new(1,Array.new(2,Array.new(3)))
=> [[[nil, nil, nil], [nil, nil, nil]]]
irb(main):002:0> a.transpose()
=> [[[nil, nil, nil]], [[nil, nil, nil]]]
irb(main):003:0>

こちらは、転換されたようです。一次元と二次元が入れ替わりました。
さらに四次元配列で確認してみます。

irb(main):001:0> a=Array.new(1,Array.new(2,Array.new(3,Array.new(4))))
=> [[[[nil, nil, nil, nil], [nil, nil, nil, nil], [nil, nil, nil, nil]], [[nil, nil, nil, nil], [nil, nil, nil, nil], [nil, nil, nil, nil]]]]
irb(main):002:0> a.transpose()
=> [[[[nil, nil, nil, nil], [nil, nil, nil, nil], [nil, nil, nil, nil]]], [[[nil, nil, nil, nil], [nil, nil, nil, nil], [nil, nil, nil, nil]]]]
irb(main):003:0>

転換されるのは一次元二次元だけで、三次元以降は変化ないようです。

Fortran90 complex


複素数を使います。
\[
\sqrt{ -2 }=
\]
この式を虚数を使って表すと、
\[
\sqrt{ 2 }i
\]
ですね。

複素数をFortranで使うには、complexを使います。

complex 変数名/(実数, 虚数)/

と定義します。

実数-2の平方根を求めるプログラムです。

takk@deb9:~$ cat func-complex-1.f90
      program main
        complex a/(-2.0, 0)/

        write(*,*) a
        write(*,*) sqrt(a)

        stop
      end
takk@deb9:~$

実行すると、実数部が0で虚数部が1.414(√2)…となりました。

takk@deb9:~$ gfortran func-complex-1.f90
takk@deb9:~$ ./a.out
 ( -2.00000000    ,  0.00000000    )
 (  0.00000000    ,  1.41421354    )
takk@deb9:~$

虚数-2の2乗を求めるプログラムです。

takk@deb9:~$ cat func-complex-2.f90
      program main
        complex a/(0, 2.0)/

        write(*,*) a
        write(*,*) a * a

        stop
      end
takk@deb9:~$ vi func-complex-2.f90
takk@deb9:~$ gfortran func-complex-2.f90
takk@deb9:~$ ./a.out
 (  0.00000000    ,  2.00000000    )
 ( -4.00000000    ,  0.00000000    )

takk@deb9:~$

実数が-4、虚数0となりました。
定数から複素数計算することもできます。

takk@deb9:~$ cat func-complex-3.f90
      program main

        write(*,*) sqrt((-4.0, 0))
        write(*,*) (0, 2)**2

        stop
      end
takk@deb9:~$ gfortran func-complex-3.f90
takk@deb9:~$ ./a.out
 (  0.00000000    ,  2.00000000    )
 ( -4.00000000    ,  0.00000000    )
takk@deb9:~$

Pythonの結果とも比べてみます。

takk@deb9:~$ python
Python 2.7.13 (default, Jan 19 2017, 14:48:08)
[GCC 6.3.0 20170118] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import cmath
>>> cmath.sqrt(-4)
2j
>>> 2j ** 2
(-4+0j)
>>>

Pythonを使うと随分計算が分かりやすくなったものだなあと思います。

コメント

タイトルとURLをコピーしました