! Program Example from "Essential Fortran 90 & 95" by Loren P. Meissner
! Copyright 1996. Copying for sale requires permission from the author.
! Otherwise, distribution is permitted if these three lines are included.

  module F23M
    implicit none
    public :: Mat_Exp
    integer, parameter, public :: P = 4, K = selected_real_kind( P )
  contains
    function Mat_Exp( M ) result( Mat_Exp_R )
      real (kind = K), dimension(:, :), intent(in) :: M
      real (kind = K), dimension(size( M, dim = 1 ), size( M, dim = 1 )) :: Mat_Exp_R
      integer :: I, N, D
      real (kind = K) :: T_Sign
      real (kind = K), allocatable, dimension(:, :) :: Term
    ! start function Mat_Exp
      Mat_Exp_R = 0.0_K      ! Initialize power series.
      T_Sign = 1.0
      D = size( M, dim = 1 )
      if (sum (M ** 2) < 0.5_K) then
        allocate( Term(D, D) )
        Term = 0.0_K  ! Initialize Term to identity matrix.
        do I = 1, D
          Term(I, I) = 1.0_K
        end do
        do N = 0, 100  ! Return when Term is small or exit after 100 terms.
          Mat_Exp_R = Mat_Exp_R + Term  ! Exp( M ) power series
          if (sum( Term ** 2 ) < 2.0_K * spacing( 1.0_K )) then
            return  ! Converged OK.
          else
            Term = T_Sign * matmul( Term, M ) / real( N + 1, kind = K )  ! Calculate next Term.
            T_Sign = - T_Sign
          end if
        end do
        write (unit = *, fmt = *) " Slow convergence. "
      end if
      return
    end function Mat_Exp
  end module F23M

