Matrix-exponentiation operator: Difference between revisions

=={{header|Fortran}}==

{{works with|Fortran|90 and later}}

<lang fortran>module matmod

implicit none

! Overloading the ** operator does not work because the compiler cannot

! differentiate between matrix exponentiation and the elementwise raising

! of an array to a power therefore we define a new operator

interface operator (.matpow.)

module procedure matrix_exp

end interface

contains

function matrix_exp(m, n) result (res)

real, intent(in) :: m(:,:)

integer, intent(in) :: n

real :: res(size(m,1),size(m,2))

integer :: i

if(n == 0) then

res = 0

do i = 1, size(m,1)

res(i,i) = 1

end do

return

end if

res = m

do i = 2, n

res = matmul(res, m)

end do

end function matrix_exp

end module matmod

program Matrix_exponentiation

use matmod

implicit none

integer, parameter :: n = 3

real, dimension(n,n) :: m1, m2

integer :: i, j

m1 = reshape((/ (i, i = 1, n*n) /), (/ n, n /), order = (/ 2, 1 /))

do i = 0, 4

m2 = m1 .matpow. i

do j = 1, size(m2,1)

write(*,*) m2(j,:)

end do

write(*,*)

end do

end program Matrix_exponentiation</lang>

Output

<pre> 1.00000 0.00000 0.00000

0.00000 1.00000 0.00000

0.00000 0.00000 1.00000

1.00000 2.00000 3.00000

4.00000 5.00000 6.00000

7.00000 8.00000 9.00000

30.0000 36.0000 42.0000

66.0000 81.0000 96.0000

102.000 126.000 150.000

468.000 576.000 684.000

1062.00 1305.00 1548.00

1656.00 2034.00 2412.00

7560.00 9288.00 11016.0

17118.0 21033.0 24948.0

26676.0 32778.0 38880.0</pre>