2つのベクトル[1,2,3]と[4,5,6]から外積(成分ごとの掛け算)により、3x3行列を作る例題を考えてみる。spreadを利用するなら
fortranなら
program test
real(8),allocatable:: aaa(:,:),bbb(:,:),ccc(:,:)
allocate(aaa,source=spread([1d0,2d0,3d0],dim=2,ncopies=1))
allocate(bbb,source=spread([4d0,5d0,6d0],dim=2,ncopies=1))
allocate(ccc,source=matmul(aaa,transpose(bbb)))
write(*,*)'aaa=',aaa
write(*,*)'bbb=',bbb
write(*,*)'outer product: aaa*bbb',shape(ccc)
write(*,*) ccc
end program test
とすればインデックス操作なしで作れる。doなどを使うとdoのレンジ指定、cccのアロケーションサイズ指定などで間違いそうな点が何点か生じる。spreadというあまりポピュラーでない関数を使わないと次元を増やせない点は面倒(このような場合、a(3)をa(3,1)にしたいだけ。unsqueezeしたいだけ。spreadって変にオーバースペック)。
とにかくfortranコードにおいても可能な限りインデックス操作やdo loopの数は減らしていくべき(無理なことをして難解になるのはいけない)。階層化をより浅くしていく、極端にはpytorchのごとくdo文を極力使わない書き方が好ましい。
- allocateする時、moldなど使わずsourceを用いて直接にshapeごとallocateすると間違いが生じにくいと思う。
- 最近のfortranにおいては、配列のshape matchingは極力自動的に行わせるようにする方が良い。配列のレンジやインデックスの直接操作はできるだけ避ける。
- fortranは一次元配列、というのは過去の常識。もちろん知っていないと困るが、shapeがとれる形になっているかどうか?shapeが引き渡せてるか?なども気にしないといけない。
- そのためにも独立した関数定義は極力避けてすべての関数、サブルーチンをモジュールに所属させる方向がよい(モジュールの使い方には主に2種類ー中途半端な使い方はまずい)。自動でインターフェイス(*.mod)を作ってくれるのはとてもありがたい。
*(原則的に自動的に決定されるような)インターフェイス文は人間がかくべきでない。仕方ない場合は仕方ないが。 - ところどころにassertを入れる(昔はsanity checkと言ってたが)。
- 情報の二重記述が安全側に働くことは稀である。たとえばmodule内でallocateしてデータ生成した場合、それを読むにはuseで直接に読むこと。配列サイズもそこから取る。subroutine冒頭でのデータ授受は極力避けるようにしていったほうが良い。
- fortranで構造体や構造体を使ったClassは極力用いていない。
pytorchならbroadcastを使えるので、以下のような感じで簡単に作れる。
import torch
# 2つのベクトルを定義
vector_a = torch.tensor([1, 2, 3])
vector_b = torch.tensor([4, 5, 6])
# ベクトルの次元を調整し、broadcastを利用
matrix = vector_a.unsqueeze(1) * vector_b
print(matrix)
これを考えると、fortranならインプリシットなbroadcastは避けることにして、以下のように書けたらいいのにな、とは思います。
(こう書けたらいいのにな)
program test
real(8),allocatable:: aaa(:,:),bbb(:,:),ccc(:,:)
allocate(ccc,source=matmul([1d0,2d0,3d0].unsqueeze(-1),[4d0,5d0,6d0].unsqueeze(1)))
write(*,*)'outer product: ',shape(ccc)
write(*,*) ccc
end program test