Working implementation of expm.

Author: loiseaujc

src/CMakeLists.txt

@@ -32,16 +32,17 @@ set(fppFiles
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp
     stdlib_linalg_eigenvalues.fypp
-    stdlib_linalg_solve.fypp    
+    stdlib_linalg_solve.fypp
     stdlib_linalg_determinant.fypp
     stdlib_linalg_qr.fypp
     stdlib_linalg_inverse.fypp
     stdlib_linalg_pinv.fypp
     stdlib_linalg_norms.fypp
     stdlib_linalg_state.fypp
-    stdlib_linalg_svd.fypp 
+    stdlib_linalg_svd.fypp
     stdlib_linalg_cholesky.fypp
     stdlib_linalg_schur.fypp
+    stdlib_linalg_matrix_functions.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp

src/stdlib_linalg.fypp

@@ -28,6 +28,7 @@ module stdlib_linalg
   public :: eigh
   public :: eigvals
   public :: eigvalsh
+  public :: expm
   public :: eye
   public :: inv
   public :: invert
@@ -1678,6 +1679,19 @@ module stdlib_linalg
       #:endfor            
   end interface mnorm
 
+!> Matrix exponential: function interface
+  interface expm
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module function expm_${ri}$(A, order) result(E)
+        ${rt}$, intent(in) :: A(:, :)
+        !! On entry, the original matrix. On exit, its exponential.
+        integer(ilp), optional, intent(in) :: order
+        !! Order of the rational approximation.
+        ${rt}$, allocatable :: E(:, :)
+    end function expm_${ri}$
+    #:endfor
+  end interface expm
+
 contains
 
 

src/stdlib_linalg_matrix_functions.fypp

@@ -0,0 +1,68 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_matrix_functions
+     use stdlib_linalg_constants
+     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
+     implicit none
+
+    contains
+
+    #:for rk,rt,ri in RC_KINDS_TYPES 
+    module function expm_${ri}$(A, order) result(E)
+        ${rt}$, intent(in) :: A(:, :)
+        integer(ilp), optional, intent(in) :: order
+        ${rt}$, allocatable :: E(:, :)
+
+        ${rt}$, allocatable :: A2(:, :), Q(:, :), X(:, :)
+        real(${rk}$) :: a_norm, c
+        integer(ilp) :: n, ee, k, s
+        logical(lk) :: p
+        integer(ilp) :: p_order
+
+        ! Deal with optional args.
+        p_order = 10 ; if (present(order)) p_order = order
+
+        n = size(A, 1)
+
+        ! Compute the L-infinity norm.
+        a_norm = mnorm(A, "inf")
+
+        ! Determine scaling factor for the matrix.
+        ee = int(log(a_norm) / log(2.0_${rk}$)) + 1
+        s = max(0, ee+1)
+
+        ! Scale the input matrix & initialize polynomial.
+        A2 = A / 2.0_${rk}$**s
+        X = A2
+
+        ! Initialize P & Q and add first step.
+        c = 0.5_${rk}$
+        E = eye(n, mold=1.0_${rk}$) ; E = E + c*A2
+
+        Q = eye(n, mold=1.0_${rk}$) ; Q = Q - c*A2
+
+        ! Iteratively compute the Pade approximation.
+        p = .true.
+        do k = 2, p_order
+            c = c*(p_order - k + 1) / (k * (2*p_order - k + 1))
+            X = matmul(A2, X)
+            E = E + c*X
+            if (p) then
+                Q = Q + c*X
+            else
+                Q = Q - c*X
+            endif
+            p = .not. p
+        enddo
+
+        E = matmul(inv(Q), E)
+        do k = 1, s
+            E = matmul(E, E)
+        enddo
+
+        return
+    end function
+    #:endfor
+
+end submodule stdlib_linalg_matrix_functions

test/linalg/CMakeLists.txt

@@ -16,6 +16,7 @@ set(
   "test_linalg_svd.fypp"
   "test_linalg_matrix_property_checks.fypp"
   "test_linalg_sparse.fypp"
+  "test_linalg_expm.fypp"
 )
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
 
@@ -35,3 +36,4 @@ ADDTEST(linalg_schur)
 ADDTEST(linalg_svd)
 ADDTEST(blas_lapack)
 ADDTEST(linalg_sparse)
+ADDTEST(linalg_expm)

test/linalg/test_linalg_expm.fypp

@@ -0,0 +1,90 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+! Test Schur decomposition
+module test_linalg_expm
+    use testdrive, only: error_type, check, new_unittest, unittest_type
+    use stdlib_linalg_constants
+    use stdlib_linalg, only: expm, eye, mnorm
+
+    implicit none (type,external)
+    
+    public :: test_expm_computation
+
+    contains
+
+    !> schur decomposition tests
+    subroutine test_expm_computation(tests)
+        !> Collection of tests
+        type(unittest_type), allocatable, intent(out) :: tests(:)
+        
+        allocate(tests(0))
+        
+        #:for rk,rt,ri in RC_KINDS_TYPES
+        tests = [tests, new_unittest("expm_${ri}$",test_expm_${ri}$)]
+        #:endfor
+
+    end subroutine test_expm_computation
+
+    !> Matrix exponential with analytic expression.
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    subroutine test_expm_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ! Problem dimension.
+        integer(ilp), parameter :: n = 5, m = 6
+        ! Test matrix.
+        ${rt}$ :: A(n, n), E(n, n), Eref(n, n)
+        real(${rk}$) :: err
+        integer(ilp) :: i, j
+
+        ! Initialize matrix.
+        A = 0.0_${rk}$
+        do i = 1, n-1
+            A(i, i+1) = m*1.0_${rk}$
+        enddo
+
+        ! Reference with analytical exponential.
+        Eref = eye(n, mold=1.0_${rk}$)
+        do i = 1, n-1
+            do j = 1, n-i
+                Eref(i, i+j) = Eref(i, i+j-1)*m/j
+            enddo
+        enddo
+
+        ! Compute matrix exponential.
+        E = expm(A)
+
+        ! Check result.
+        err = mnorm(Eref - E, "inf")
+        call check(error, err < (n**2)*epsilon(1.0_${rk}$), "Analytical matrix exponential.")
+        if (allocated(error)) return
+        return        
+    end subroutine test_expm_${ri}$    
+    #:endfor    
+
+end module test_linalg_expm
+
+program test_expm
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_linalg_expm, only : test_expm_computation
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [ &
+        new_testsuite("linalg_expm", test_expm_computation) &
+        ]
+
+    do is = 1, size(testsuites)
+        write(error_unit, fmt) "Testing:", testsuites(is)%name
+        call run_testsuite(testsuites(is)%collect, error_unit, stat)
+    end do
+
+    if (stat > 0) then
+        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        error stop
+    end if
+end program test_expm