t-nissieの日記： 【電脳】FortranでAPL風の直積 (outer product) またはディアド (dyadic product)

おもてのAPLのストーリーで、学部学生だったころ情報処理教育センターの
3270端末のキーボードにはAPL用の記号が書いてあったのを懐かしく思い
出しました。OSはたしかVM/CMSとかいうやつでした。すぐにUNIXに
変更になりましたが。

WikipediaのAPLを読んで、直積とかが１行で書けるのに感動しました。
こりゃあ頭のよい金融屋が思いついたことすぐAPLで書いてがっつり儲ける
わけだと思いました。

で、Fortranで直積を書いてみました。さいきんのFortranはモジュールに
するかinterface文を使えばサブルーチンに配列の長さを渡せるのです。
３×３のディアドはたまに使いますが、いままでは引数の長さが固定の
関数かサブルーチンで書いてました。

でもこれ、uとvがすごく大きかったらスタックがあふれるよね。

$ make -k && ./use_outer_product
gfortran -ffree-form -Wall -O3 -g  -c -o outer_product_module.o outer_product_module.f
gfortran -ffree-form -Wall -O3 -g  -c -o use_outer_product.o use_outer_product.f
gfortran -ffree-form -Wall -O3 -g  -o use_outer_product use_outer_product.o outer_product_module.o
  5.0  6.0  7.0
  10.0 12.0 14.0
  15.0 18.0 21.0
  20.0 24.0 28.0

Makefile:

# -*-Makefile-*-
# Time-stamp: <2009-06-12 10:51:57 t-nissie>
# Author: t-nissie
##
FC=gfortran
FFLAGS=-ffree-form -Wall -O3 -g
 
use_outer_product: use_outer_product.o outer_product_module.o
    $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
 
use_outer_product.o: outer_product_module.o
 
clean:
    rm -f core *.mod *.o use_outer_product a.out

outer_product_module.f:

! outer_product_module.f -*-f90-*-
! Time-stamp: <2009-06-12 10:40:47 t-nissie>
! Author: t-nissie
!!
module outer_product_module
implicit none
contains
  function outer_product(u,v)
    implicit none
    real*8, intent(in) :: u(:)
    real*8, intent(in) :: v(:)
    real*8             :: outer_product(size(u),size(v))
    integer m,n,i,j
    m=size(u)
    n=size(v)
    do j = 1, n
       do i = 1, m
          outer_product(i,j) = u(i)*v(j)
       end do
    end do
  end function outer_product
end module outer_product_module

use_outer_product.f:

! use_outer_product.f -*-f90-*-
! Time-stamp: <2009-06-12 10:35:59 t-nissie>
! Author: t-nissie
!!
program use_outer_product
  use outer_product_module
  implicit none
  real*8 u(4)
  real*8 v(3)
  real*8 o(4,3)
  integer i,j
  data u/1.0d0, 2.0d0, 3.0d0, 4.0d0/
  data v/5.0d0, 6.0d0, 7.0d0/
  o = outer_product(u,v)
  do i=1,4
     write(6,'(3f5.1)') (o(i,j),j=1,3)
  end do
end program use_outer_product

参考にしたのはここ。