FORTRAN77を学習してプログラミングできるようになった気がしますので、Fortran90の方へ移っていきたいと思います。
Contents
- Fortran90 ビルド、継続行
- Fortran90 do
- Fortran90 do exit
- Fortran90 do cycle
- Fortran90 大文字小文字区別
- Fortran90 定数
- Fortran90 data
- Fortran90 type
- Fortran90 where
- Fortran90 where その2
- Fortran90 where その3
- Fortran90 多次元配列
- Fortran90 多次元配列 その2
- Fortran90 多次元配列 その3
- Fortran90 行列
- Fortran90 行列の積
- Fortran90 sum
- Fortran90 sum その2
- Fortran90 product
- Fortran90 shape,maxval,minval,reshape
- Fortran90 cshift,eoshift
- Fortran90 merge
- Fortran90 spread
- Fortran90 transpose
- Fortran90 complex
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を使うと随分計算が分かりやすくなったものだなあと思います。
コメント