2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

fortran spreadの使い方:2つのベクトルの外積で行列を作る

Last updated at Posted at 2025-02-09

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
2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?