線形代数計算

行列やベクトルの定義

Fortranでは行列やベクトルの定義の方法についてです。配列の定義方法について比較してきましょう。Fortranでは配列の定義は

real(8)::d(1:3,1:4)

のような形で書きますね（他の書き方は省略）。この場合はdは倍精度実数が要素の2次元配列です。この配列の中身は決まっていませんので、通常は

d = 0d0

などと初期化します。Juliaの場合、0で初期化した2次元配列が必要な場合

d = zeros(Float64,3,4)

と書きます。もし、Fortranと同じように初期化されていない配列が欲しい場合には、

d = Array{Float64,2}(undef,3,4)

と書きます。Array{Float64,2}の2は配列の次元を示し、(undef,3,4)の3,4はそれぞれの次元の要素数を示します。

倍精度複素数の3次元配列が欲しい場合は、Fortranであれば

complex(8)::f(1:2,1:3,1:5)

で定義できます。Juliaであれば

f = zeros(ComplexF64,2,3,5)

となります。

Fortranではallocatableを使えば

real(8),allocatable::d(:,:)
allocate(d(1:3,1:4))

というような形で、実行時に動的に配列の形状を決めることができていました。Juliaのzeros(Float64,3,4)やArray{Float64,2}(undef,3,4)などはこれと同じように動的に配列の形状を決めていると考えても構いません。両者はよく似ています。

Fortranの場合、配列の始まりと終わりを好きなように決めることができました。例えば、

real(8)::d(1:2,3:9)

というdを定義すると、1次元目の添字の範囲は1から2、2次元目の添字の範囲が3から9の配列を定義することができます。Juliaではデフォルトではこのような添字の範囲を設定することはできません。しかし、JuliaのパッケージにはOffsetArraysというパッケージがありまして、これを使うことで同じことが可能となります。例えば、上のような配列の場合には

using OffsetArrays
c = Array{Float64,2}(undef,2,7)
d = OffsetArray(c, 1:2, 3:9)

とすればdの添字の範囲はFortranと同じになります。なお、OffsetArraysを使って通常の配列cからdを作成する時にメモリのコピーは発生しておりませんので、追加のコストはほとんどかかりません。

初期値について

Juliaでは、作った配列の初期値をあらかじめ0とする

a = zeros(Float64,2,3)

の他に、全部1にする

b = zeros(Int64,2,3)

があります。また、乱数での初期化として

c = rand(2,3)

も可能です。乱数は0から1の間の数が出ます。なお、乱数のシードは

using Random
n = 123
Random.seed!(n)

で固定できます。nを変えると異なる乱数列が出ます。nが同じの場合には常に同じ数列が出力されますから、デバッグなどでは重要です。

基本的な演算

Fortranも行列やベクトルの演算に強いですが、Juliaも強いです。例えば、Fortranでは

y = matmul(A,x)
C = matmul(A,B)

などで、$y = Ax$や$C = AB$が計算できます。ここで大文字は行列、小文字はベクトルを想定しました。Juliaでは

y = A*x
C = A*B

とただの掛け算記号でかけます。もちろん、足し算もできて、ベクトル同士の足し算は

y = A*x + b

でできます。

線形代数関連

Fortranで連立方程式を解いたり逆行列を計算したり特異値分解をしたり固有値を求めたりする場合には、LAPACKをインストールして使うと思います。Juliaでは標準機能でこれらを実行できます。まず、線形代数関連のパッケージを使うため、

using LinearAlgebra

とします。連立方程式$b = Ax$の$x$を求めたい場合には

y = rand(2,2)
b = rand(2,2)
x = A \ b

ですぐ求まります。