fortranを用いた行列の掛け算について

行列の掛け算についてなのですが、ある行列aを２乗した行列bを求める場合は以下のようなプログラムを書けば出来たのですが、これを３乗以上に拡張するためにはどうしたらよいのでしょうか？

______do i=1,3
________do j=1,3
__________b(i,j)=0.D0
____________do k=1,3
______________b(i,j)=b(i,j)+a(i,k)*a(k,j)
____________enddo
________enddo
______enddo

No.1 回答者： himajin100000 回答日時： 2009/10/20 07:22

program hello

implicit none
real,allocatable,dimension(:,:):: d
real,allocatable,dimension(:,:):: x

integer::i

allocate(d(3,3))
allocate(x(3,3))

d(1,1) = 5
d(1,2) = 4
d(1,3) = 4

d(2,1) = 7
d(2,2) = 2
d(2,3) = 4

d(3,1) = 7
d(3,2) = 6
d(3,3) = 2

x = mpow(d,3);

do i = 1,ubound(x,2)
print *, x(i,:)
end do

contains

function mmulti(x,y)

real,intent(in),allocatable,dimension(:,:)::x
real,intent(in),allocatable,dimension(:,:)::y
real,allocatable,dimension(:,:)::mmulti

integer::i
integer::j
integer::k

allocate(mmulti(ubound(x,1),ubound(y,2)))

do i = 1,ubound(x,1)
do j = 1,ubound(y,2)
mmulti(i,j) = 0
do k = 1,ubound(x,2)
mmulti(i,j) = mmulti(i,j) + x(i,k) * y(k,j)
end do
end do
end do

end function mmulti

function mpow(x,n)

real,intent(in),allocatable,dimension(:,:)::x
integer,intent(in)::n
real,allocatable,dimension(:,:)::mpow

integer::i

allocate(mpow(ubound(x,1),ubound(x,1)))

mpow = x
do i = 1,n - 1
mpow = mmulti(mpow,x);
end do


end function mpow

end program hello


全部書いてみた。g95でコンパイルして

1077. 692. 620.
1085. 684. 620.
1211. 804. 684.

の結果を得て、Excelとの一致も確認しました。

表示用のサブルーチンに分けようとしてうまくいかなかったり、
n乗で与えられるnが負でないときや
与えられた二つの行列x,yについて,xの行数とyの列数が等しくないときや
xが正方行列でないときに
エラー出さないといけないんですが，
方法をすぐに調べられなかったので断念しました。