行列(2次元配列)

次のプログラムは,以下に示す2行3列の行列 A, B の和を求め,その結果を表示するプログラムです.

A = ( 1 2 3 )         B = ( 6 3 4 )
    ( 4 5 6 )             ( 5 1 2 )
program work0701
  implicit none
  integer :: i, j
  integer :: ma(2, 3) = reshape( (/ 1, 4, 2, 5, 3, 6 /), (/ 2, 3 /) )   ! 行列 A
  integer :: mb(2, 3) = reshape( (/ 6, 5, 3, 1, 4, 2 /), (/ 2, 3 /) )   ! 行列 B
  integer :: mc(2, 3) = 0                                               ! 行列 C

  do i = 1, 2
     do j = 1, 3
        mc(i, j) = ma(i, j) + mb(i, j)
     end do
  end do

  do i = 1, 2
     write(*,*) mc(i, 1:3)
  end do

  stop
end program work0701

配列 ma, mb, mc は「2次元配列」と呼ばれます. 要素数が 2×3 である2次元配列 ma(2, 3) は,次のように integer 型の変数が並んだものと考えればよいでしょう.

┌────┬────┬────┐
│ma(1, 1)│ma(1, 2)│ma(1, 3)│
├────┼────┼────┤
│ma(2, 1)│ma(2, 2)│ma(2, 3)│
└────┴────┴────┘

このように,2次元配列は縦×横の2次元の表とみなすことができ,各要素を指定するには2つの添字が必要になります.

2次元以上の配列を「多次元配列」といいます. Fortran 90/95 では7次元配列まで扱うことができます.

配列要素順序

概念的には2次元の表である2次元配列も,コンピュータの内部では次のような順番で要素が連続して並んでいます.

┌────┬────┬────┬────┬────┬────┐
│ma(1, 1)│ma(2, 1)│ma(1, 2)│ma(2, 2)│ma(1, 3)│ma(2, 3)│
└────┴────┴────┴────┴────┴────┘

この順序のことを「配列要素順序」といいます. Fortran では左の添字ほど速く変化するような順番になっています. これを column-major order といいます.

2次元配列の初期化

2次元配列を初期化するには reshape 関数を利用します. reshape 関数の第1引数は初期化データを並べた1次元配列,第2引数は配列の「形状」(2行3列の行列なら (/ 2, 3 /) など.後述)です. 第1引数の各要素は配列要素順序にしたがって値を並べることに注意しましょう.

練習問題

以下に示す行列 A, B の積 C を求めるプログラムを作成せよ.

A = ( 1 2 3 )     B = ( 1 4 )
    ( 4 5 6 )         ( 2 5 )
                      ( 3 6 )

行列 C の i, j 成分 Cij は次のように求められる.

Cij = Σ Aik Bkj      (Σ は k=1,2,3 についての和を取る)

これを各 i, j について行えばよい.

クリックして表示

添字の下限と上限

問題によっては配列の添字を 1 以外の数から始めたいことがあります. Fortran では添字の上限だけでなく下限を指定することができます. 次のように宣言すると,nx は要素数 16 の配列 nx(-10), nx(-9), …, nx(5) となります.

integer :: nx(-10:5)

下限を省略したときは 1 が指定されたものとみなされます.

2次元以上の配列でも同様です.

integer :: ma(0:1, -1:1)

と宣言すると,2次元配列 ma の要素数は 6 となります.

┌────┬────┬────┐
│ma(0,-1)│ma(0, 0)│ma(0, 1)│
├────┼────┼────┤
│ma(1,-1)│ma(1, 0)│ma(1, 1)│
└────┴────┴────┘

階数・寸法・形状

配列の次元数のことを階数 (rank),各次元の要素数のことを寸法 (extent) といいます. 次のように宣言された3次元配列では,階数が 3,各次元の寸法が 11, 6, 10 となります.

real(8) :: grid(-5:5, 0:5, 10)

また,寸法を順に並べたものを形状 (shape) といいます. 配列 grid の形状は (11, 6, 10) となります.

練習問題

次の配列の形状と要素数を求めよ.

integer :: ma(3, 3)
integer :: mb(-3:3, 1:5)
real(8) :: ta(0:9, 0:9, -1:1)

部分配列

Fortran 90 以降では,配列の要素をある範囲で参照することができます. これを「部分配列」といいます. 部分配列では,各次元の範囲を次のように指定します.

[lower]:[upper][:stride]

lower または upper が省略されると,配列宣言に用いられた下限と上限が指定されたものとみなされます. stride は数個おきに要素を参照する場合に用いられます. stride が省略された場合は 1 が仮定され,要素が連続的に参照されます. stride には負の整数を指定することもできます.

a(1:5)          ! a(1), a(2), a(3), a(4), a(5) を要素とする1次元配列
a(1:10:2)       ! a(1), a(3), a(5), a(7), a(9) を要素とする1次元配列
a(:)            ! 1次元配列 a の全体
g(1, 1:3)       ! g(1,1), g(1,2), g(1,3) を要素とする1次元配列
h(0:1, 2:3)     ! h(0,2), h(0,3), h(1,2), ..., h(1,3) を要素とする2次元配列

次のプログラムでは,1次元配列 nx のいろいろな範囲を表示しています.

program work0703
  implicit none
  integer :: i
  integer :: nx(1:9) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9 /)

  write(*,*) nx(1:3)         ! nx(1), nx(2), nx(3)
  write(*,*) nx(7:)          ! nx(7), nx(8), nx(9)
  write(*,*) nx(9:1:-2)      ! nx(9), nx(7), nx(5), nx(3), nx(1)

  stop
end program work0703

上限と下限をともに省略し,配列全体を nx(:) と書くこともできます.

練習問題

work0703 を書き換えて,配列 nx の全部の要素を逆順に表示せよ. 部分配列を利用すること.

クリックして表示

配列の基本演算

配列演算

Fortran 90 以降では加減乗除やべき乗,多くの組込み関数が配列に対応しています. ただし,配列に関する計算では

ことに注意が必要です.

配列が次のように宣言されているとします.

real(8) :: a(10, 10), b(10, 10)   ! 10x10 行列
real(8) :: x(5), y(5)             ! 5 次元ベクトル

配列全体やその一部(部分配列)を対象とした計算が可能です.

- a                 ! a の全要素に -1 をかける: - a(i, j)
a * b               ! a と b の対応する要素をかける: a(i, j) * b(i, j)
x + 1               ! x の全要素に 1 を加える: x(i) + 1
a**2                ! a の全要素を2乗する: a(i, j)**2
1 / x + a(1, 2:6)   ! 1/x(i) + a(1, i+1) (i=1, 2, ..., 5)
a == b              ! a(i, j) == b(i, j) がすべての i, j について成り立てば真,
                    ! そうでなければ偽
abs(x)              ! abs(x(i))
sin(a) * cos(b)     ! sin(a(i, j)) * cos(b(i, j))

配列に対する乗算 a * b は行列積ではないことに注意しましょう. また 1 / xa * b と書くと x や a, b がスカラーなのか配列なのか,これを見ただけでは判断できません. 配列であることを強調するために 1 / x(:), a(:,:) * b(:,:) と書く流儀もあります.

配列代入

配列演算の結果を同じ形状の配列に代入することができます.

a = 0               ! a の全要素に 0 を代入する: a(i, j) = 0
a = a + 1           ! a の全要素に 1 を加える: a(i, j) = a(i, j) + 1
a(5, 6:10) = x      ! a(5, 6:10) に x(1:5) をコピーする: a(5, i+5) = x(i)
y = sin(x)          ! y(i) = sin(x(i))

添字の番号が異なる場合,配列要素の対応関係に注意しましょう.

配列要素の対応

配列代入が行われる場合,配列要素の対応は添字の番号ではなく各次元での位置(何番目の要素であるか)で行われます. 例として,

real(8) :: a(10, 5), b(0:9, 0:4), c(11:21, 11:15)

という3つの配列を考えましょう. これらの配列は添字の範囲が異なっていますが,すべて同じ形状になっており,配列演算・配列代入が可能です.

c = a + b

という配列演算では次のように要素が対応します.

do j = 1, 5
   do i = 1, 5
      c(i+10, j+10) = a(i, j) + b(i-1, j-1)
   end do
end do

ただし,実際の計算順序(どの要素から計算が計算が行われるか)は不定です.

配列代入の注意点

配列代入は「右辺がすべて評価されてから行われる」ことに注意が必要です.

integer :: nx(5) = (/ 1, 2, 3, 4, 5 /)
nx(2:5) = nx(1:4)

を実行すると配列 nx は (/ 1, 1, 2, 3, 4 /) となります. この配列代入文を

do i = 2, 5
   nx(i) = nx(i-1)
end do

と理解してはいけません. do ループの場合,代入は逐次的に行われるため,

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

が上から順に評価され,最終的に nx は (/ 1, 1, 1, 1, 1 /) となります.

練習問題

配列 nx, ny を次のように宣言する.

integer :: nx(5) = (/ 1, 2, 3, 4, 5 /)
integer :: ny(5)

配列代入を用いて,nx の各要素を逆順に ny にコピーするプログラムを作成せよ.

クリックして表示

組込み関数

最大値と最小値 maxval, minval

Fortran 90 以降では,組込み関数 maxval(配列名), minval(配列名) を使って,配列の最大値・最小値を求めることができます. 配列は何次元配列でも構いません.

tmp = maxval(nx)    ! nx(1:5) の最大値

とすると,tmp には nx の最大値が代入されます.

練習問題

キーボードから5つの整数を入力し,その最大値を求めるプログラムを作成せよ. 組込み関数 maxval を用いること.

クリックして表示

総和 sum

組込み関数 sum を使うと配列の要素の和を求めることができます. 配列は何次元配列でも構いません.

次の例は,キーボードから5つの整数を入力し,その総和を求めるプログラムです.

program work0707
  implicit none
  integer :: nx(5)

  write(*,*) 'Input 5 integers:'
  read(*,*) nx(1:5)

  write(*,*) sum(nx)

  stop
end program work0707

sum の引数に部分配列を与えると,配列の一部の和を求めることができます. sum(nx(1:3)) とすれば,nx の1番目から3番目の要素の和となります.

練習問題

work0707 を書き換え,平均値を求めるプログラムを作成せよ.

クリックして表示

本日の課題

内積の計算

2つの3次元ベクトル x, y をキーボードから読み込み,内積 x.y を計算するプログラムを作成せよ. x, y は real(8) 型の配列とする. 配列同士の乗算 x * y と,組込み関数 sum を利用せよ.

クリックして表示

パスカルの三角形

次の下三角行列はパスカルの三角形と呼ばれる.

1   0   0   0    …     0
1   1   0   :           :
1   2   1   0
1   3   3   1   0
1   4   6   4   1   0   :
1   5  10  10   5   1   0
1   6  15  20  15   6   1

この行列を配列 p(1:9, 1:9) に作れ.

ヒント: 上三角行列はすべて 0,第 1 列をすべて 1 とすると,他の要素はすべて真上と左上の要素の和になる.

クリックして表示

Newton 法

キーボードから実数 a を与え,方程式 tan(x) = a をNewton 法で解くプログラムを作成せよ. 漸化式は次で与えられる.

Xn+1 = Xn - (tan(Xn) - a) * cos(Xn)**2

初期値 X0 は 0 とせよ.

クリックして表示