能書き

ここでは、配列の高度な操作を学びます。
ここに書いてあること以外にも様々なルールがありますが、 その中から「使える」ものを厳選しました。

program hairetsu
 !------------------------
 ! Fortran template
 !------------------------
 use bacs
 implicit none


end program

定義：宣言

後から大きさを決められる配列を宣言する

 型, allocatable :: ベクトル名(:), 行列名(:,:)

本文中で配列の大きさを決める

 allocate(配列名(index1の上限[,index2の上限,...]), [配列名2...])

学習

「行列の大きさ」は、物理の問題では「系の大きさ」に相当する場合がしばしばあります。
系の大きさを変えるたびにプログラム中の数十カ所に散らばった値を１つ１つ書き直す、 なんてことをしていては卒研がいつまでも終わりません。

たとえば大きさ『n』の系を扱うときは次のようにします。

 integer :: n
 real(kind=double), allocatable :: vector(:), matrix(:,:)

 read *, n
 allocate( vector(n), matrix(n,n) )

たとえばnに『4』を入れ『allocate』することで、大きさ『4』のベクトル、行列として使えるようになります。 大きさを２倍にするには、実行時に『8』と入れれば良いだけです。
このように「柔軟」な書き方をしておくのが、上手なプログラミングのポイントです。

学習： 配列代入

 配列名 = スカラー変数名

『配列名』のすべての要素に『スカラー変数名』を代入する。

 配列名1 = 配列名2

『配列名1』の各要素に『配列名2』の対応する要素を代入する。
配列同士の形状が等しくなければならない。

学習

「Amatの要素をすべて1にする」「AmatをBmatにコピーする」
は、普通に書くと

 do i=1,n
    do j=1,n
       Amat(i,j) = 1.0d0
    end do
 end do

 do i=1,n
    do j=1,n
       Bmat(i,j) = Amat(i,j)
    end do
 end do

ですが、配列代入文を使うと

 Amat = 1.0d0

 Bmat = Amat

と書けます。こっちの方が簡単キレイですね。

定義： 四則演算・初等関数の配列版

 配列名1 二項演算子 配列名2

『配列名1』と『配列名2』の対応する各要素同士を演算する。
ただし、配列同士の形状が等しくなければならない。

 print *, 配列名

『配列名』の全要素を横に並べて表示する。

 初等関数名(配列名)

『配列名』の各要素に対する『初等関数』を要素とする配列を返す。

学習

行列の和、差は

 Cmat = Amat + Bmat

 Cmat = Amat - Bmat

と書けて便利です。
積と商は、まず使わないでしょう。むしろ、

 『Amat*Bmat』は行列の積ではない
 『1/Amat』『Amat**(-1)』は逆行列ではない

ということを理解しておいてください。

 print *, Amat

とすると、『Amat』の中身がズラズラと表示されます。 上手く表示させるには、工夫が必要です。
配列の初等関数は

 log(Amat)

くらいは使い道があるかもしれませんが・・・これもむしろ、

 『exp(Amat)』はeの肩に行列を乗せたものではない

ということを理解しておいてください。

実践

単位行列を定義するプログラムを作ってみましょう。

 integer :: i, j, n
 real(kind=double), allocatable :: Imat(:,:)
 
 print *, "n = ?"
 read *, n
 allocate( Imat(n,n) )
 Imat = zero
 do i=1,n
    Imat(i,i) = one
 end do
 print *, "Identity matrix"
 do i=1,n
    print *, Imat(i,1:n)
 end do

基本編でも一度作りましたが（参照）、 今回の方が洗練され、実用的になっています。 さて、コピー＆ペーストした人は気づかなかったかもしれませんが、

    print *, Imat(i,1:n)

という部分があります。これについては次で。

定義： 部分配列

 配列名(下限:上限 [,...])

配列の一部（『下限』から『上限』まで）からなる配列を使う。
『上(下)限』は省略すると自動的に「配列が定義された上(下)限」になる。

学習

たとえば「行列をブロックで扱う」ときに使います。 (2n)x(2n)行列の右上のnxnブロックだけ使いたい時は、次のようにします。

 A2mat(1:n,1:n) = Amat(1:n,(n+1):(2*n))

「行列からベクトルを取り出す」こともできます。
たとえば行列の「i行」を使いたいときは、次のように行のindexのみ固定します。

 Amat(i,1:n)

これをdoループを使って１行ごとにprintすれば、 行列の形で行列を表示させることができる、というわけです。

『Amat』がnxn行列の場合、下の３つは全く同じ挙動になります。

 Amat          = zero     ! パターン１
 Amat(:,:)     = zero     ! パターン２
 Amat(1:n,1:n) = zero     ! パターン３

パターン１だと配列の次元が分からないので、パターン２か３を使うのがお勧めです。

学習

「下限:上限」という書き方は、配列の宣言時にも使えます。

 allocate(Bmat(1:n,1:n))

は、次の命令と同等です。

 allocate(Bmat(n,n))

例えば以下のようにすると、「0からn+1までの要素」を持つ配列を定義できます。

 allocate(Bmat(0:(n+1),0:(n+1)))

実践

次のルーチンは何をするものか、解読して下さい。

  Cmat(:,:) = cZero   ! cZero = (0.0d0, 0.0d0)
  do j=1,m
     do k=1,n
        Cmat(1:m,j) = Cmat(1:m,j) + Amat(1:m,k)*conjg(Bmat(j,k))
     end do
  end do

理解できたら、下のような形で使ってみましょう。

 program hairetsu
 use bacs
 implicit none
 complex(kind=double), allocatable  :: Amat(:,:), Bmat(:,:), Cmat(:,:)
 integer :: m,n

  m = 2
  n = 3
  allocate(...) ! 自分で設定すること

  ! --- define matrices ---   注：ここの値は適当に決める。
  Amat(:,:)   = cOne    ! cOne = (1.0d0, 0.0d0)
  Bmat(:,:)   = cZero
  Bmat(1,1:n) = img     ! img  = (0.0d0, 1.0d0) を１行目だけに入れてみる

  Cmat = ...   ! 上のルーチンを貼り付ける
 
  ! --- print C matrix ---   ! ループを使ってうまく表示してみよう。
  print *, "Real part of C matrix"
  print *, dble(Cmat)        ! 実数部分を(倍精度で)取り出す
  print *, "Imagnary part of C matrix"
  print *, aImag(Cmat)       ! 虚数部分を取り出す

 end program

予想した結果と同じになりましたか？
配列の扱い方は最初は難しいと感じるかもしれませんが、徐々に慣れていって下さい。

紹介： 配列関数

各要素ではなく「配列そのもの」に作用する関数がいくつか用意されている。
以下に、よく使うものを列挙する。

行列の積（戻り値：行列）
 matmul(行列1,行列2)

ベクトルの内積（スカラー）
 dot_product(ベクトル1,ベクトル2)

転置行列（行列）
 transpose(行列)

複素共役行列（行列）
 conjg(行列)

配列の全成分の和（スカラー）
 sum(配列)

ベクトルの要素の中の最大値、最小値（スカラー）
 maxval(ベクトル)、 minval(ベクトル)

配列の[i番目の添字の]成分の数（スカラー）
 size(配列[,i])

まとめ

初心者は無闇に格好つけず、まず要素の指定がされた「行列の成分」の扱いをマスターしてください。
そして理解できたものから少しずつ、今回の内容を取り入れていくと良いでしょう。
ちなみに、Fortran90はC++など他の言語と比べ、配列の扱いが非常に洗練されています。
そのため、特に大規模な行列計算が求められるスパコン向きの言語と言えます。