Fortran90(行列の積)

前回は行列の内積になってしまったので、通常の積をしてみます。この式でしたね。
\[
\left(
\begin{array}{cc}
a11 & a12 \\
a21 & a22
\end{array}
\right)
\times \left(
\begin{array}{cc}
b11 & b12 \\
b21 & b22
\end{array}
\right)
= \left(
\begin{array}{cc}
a11 \times b11 + a12 \times b21 & a11 \times b12 + a12 \times b22 \\
a21 \times b11 + a22 \times b21 & a21 \times b12 + a22 \times b22
\end{array}
\right)
\]
式をそのままプログラムしてみます。

takk@deb9:~$ cat matrix-mul.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /1,2, 3,4/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)

        c(1,1) = a(1,1) * b(1,1) + a(1,2) * b(2,1)
        c(1,2) = a(1,1) * b(1,2) + a(1,2) * b(2,2)
        c(2,1) = a(2,1) * b(1,1) + a(2,2) * b(2,1)
        c(2,2) = a(2,1) * b(1,2) + a(2,2) * b(2,2)

        write(*,*) 'c=a*b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$

実行すると、

takk@deb9:~$ gfortran matrix-mul.f90
takk@deb9:~$ ./a.out
 a=
           1           3
           2           4
 b=
          10          30
          20          40
 c=a*b
          70         150
         100         220
takk@deb9:~$

行列の積をするmatmulという組み込み関数があるらしいので、使ってみます。

takk@deb9:~$ cat matrix-mul2.f90
      program main
        integer a(2,2),b(2,2),c(2,2)
        data a /1,2, 3,4/
        data b /10,20, 30,40/

        write(*,*) 'a='
        write(*,*) (a(1,i),i=1,2)
        write(*,*) (a(2,i),i=1,2)

        write(*,*) 'b='
        write(*,*) (b(1,i),i=1,2)
        write(*,*) (b(2,i),i=1,2)

        c = matmul(a,b)

        write(*,*) 'c=a*b'
        write(*,*) (c(1,i),i=1,2)
        write(*,*) (c(2,i),i=1,2)

        stop
      end
takk@deb9:~$

実行すると、

takk@deb9:~$ gfortran matrix-mul2.f90
takk@deb9:~$ ./a.out
 a=
           1           3
           2           4
 b=
          10          30
          20          40
 c=a*b
          70         150
         100         220
takk@deb9:~$

c=a*bの答えが同じになりました。

Leave a Reply

Your email address will not be published. Required fields are marked *

CAPTCHA