Matrix exponential

doc/specs/stdlib_linalg.md

@@ -1884,3 +1884,72 @@ If `err` is not present, exceptions trigger an `error stop`.
 {!example/linalg/example_mnorm.f90!}
 ```
 
+## `expm` - Computes the matrix exponential {#expm}
+
+### Status
+
+Experimental
+
+### Description
+
+Given a matrix \(A\), this function computes its matrix exponential \(E = \exp(A)\) using a Pade approximation.
+
+### Syntax
+
+`E = ` [[stdlib_linalg(module):expm(interface)]] `(a [, order])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the data. It is an `intent(in)` argument.
+
+`order` (optional): Shall be a non-negative `integer` value specifying the order of the Pade approximation. By default `order=10`. It is an `intent(in)` argument. 
+
+### Return value
+
+The returned array `E` contains the Pade approximation of \(\exp(A)\).
+
+If `A` is non-square or `order` is negative, it raises a `LINALG_VALUE_ERROR`.
+
+### Example
+
+```fortran
+{!example/linalg/example_expm.f90!}
+```
+
+## `matrix_exp` - Computes the matrix exponential {#matrix_exp}
+
+### Status
+
+Experimental
+
+### Description
+
+Given a matrix \(A\), this function computes its matrix exponential \(E = \exp(A)\) using a Pade approximation.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):matrix_exp(interface)]] `(a [, e, order, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the data. If `e` is not passed, it is an `intent(inout)` argument and is overwritten on exit by the matrix exponential. If `e` is passed, it is an `intent(in)` argument and is left unchanged.
+
+`e` (optional): Shall be a rank-2 `real` or `complex` array with the same dimensions as `a`. It is an `intent(out)` argument. On exit, it contains the matrix exponential of `a`.
+
+`order` (optional): Shall be a non-negative `integer` value specifying the order of the Pade approximation. By default `order=10`. It is an `intent(in)` argument. 
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. 
+
+### Return value
+
+The returned array `A` (in-place) or `E` (out-of-place) contains the Pade approximation of \(\exp(A)\).
+
+If `A` is non-square or `order` is negative, it raises a `LINALG_VALUE_ERROR`.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_matrix_exp.f90!}
+```
+

example/linalg/CMakeLists.txt

@@ -57,3 +57,5 @@ ADD_EXAMPLE(qr)
 ADD_EXAMPLE(qr_space)
 ADD_EXAMPLE(cholesky)
 ADD_EXAMPLE(chol)
+ADD_EXAMPLE(expm)
+ADD_EXAMPLE(matrix_exp)

example/linalg/example_expm.f90

@@ -0,0 +1,18 @@
+program example_expm
+   use stdlib_linalg, only: expm
+   implicit none
+   real :: A(3, 3), E(3, 3)
+   integer :: i
+   A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+   E = expm(A)
+
+   print *, "Matrix A :"
+   do i = 1, 3
+      print *, A(i, :)
+   end do
+
+   print *, "Matrix exponential E = exp(A):"
+   do i = 1, 3
+      print *, E(i, :)
+   end do
+end program example_expm

example/linalg/example_matrix_exp.f90

@@ -0,0 +1,20 @@
+program example_expm
+   use stdlib_linalg, only: matrix_exp
+   implicit none
+   real :: A(3, 3), E(3, 3)
+   integer :: i
+
+   print *, "Matrix A :"
+   do i = 1, 3
+      print *, A(i, :)
+   end do
+
+   A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+   call matrix_exp(A) ! In-place computation.
+   ! For out-of-place, use call matrix_exp(A, E).
+
+   print *, "Matrix exponential E = exp(A):"
+   do i = 1, 3
+      print *, E(i, :)
+   end do
+end program example_expm

src/CMakeLists.txt

@@ -46,6 +46,7 @@ set(fppFiles
     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_constants.fypp

@@ -71,6 +71,7 @@ module stdlib_constants
     #:for k, t, s in R_KINDS_TYPES
     ${t}$, parameter, public :: zero_${s}$ = 0._${k}$
     ${t}$, parameter, public :: one_${s}$  = 1._${k}$
+    ${t}$, parameter, public :: log2_${s}$ = log(2.0_${k}$)
     #:endfor
     #:for k, t, s in C_KINDS_TYPES
     ${t}$, parameter, public :: zero_${s}$ = (0._${k}$,0._${k}$)

src/stdlib_linalg.fypp

@@ -28,6 +28,7 @@ module stdlib_linalg
   public :: eigh
   public :: eigvals
   public :: eigvalsh
+  public :: expm, matrix_exp
   public :: eye
   public :: inv
   public :: invert
@@ -1678,6 +1679,107 @@ module stdlib_linalg
       #:endfor            
   end interface mnorm
 
+  !> Matrix exponential: function interface
+  interface expm
+    !! version : experimental
+    !!
+    !! Computes the exponential of a matrix using a rational Pade approximation.
+    !! ([Specification](../page/specs/stdlib_linalg.html#expm))
+    !!
+    !! ### Description
+    !!
+    !! This interface provides methods for computing the exponential of a matrix
+    !! represented as a standard Fortran rank-2 array. Supported data types include
+    !! `real` and `complex`.
+    !!
+    !! By default, the order of the Pade approximation is set to 10. It can be changed
+    !! via the `order` argument that must be non-negative.
+    !!
+    !! If the input matrix is non-square or the order of the Pade approximation is
+    !! negative, the function returns an error state.
+    !!
+    !! ### Example
+    !!
+    !! ```fortran
+    !!  real(dp) :: A(3, 3), E(3, 3)
+    !!
+    !!  A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+    !!
+    !!  ! Default Pade approximation of the matrix exponential.
+    !!  E = expm(A)
+    !!
+    !!  ! Pade approximation with specified order.
+    !!  E = expm(A, order=12)
+    !! ```
+    !!
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module function stdlib_linalg_${ri}$_expm_fun(A, order) result(E)
+        !> Input matrix a(:, :).
+        ${rt}$, intent(in) :: A(:, :)
+        !> [optional] Order of the Pade approximation (default `order=10`)
+        integer(ilp), optional, intent(in) :: order
+        !> Exponential of the input matrix E = exp(A).
+        ${rt}$, allocatable :: E(:, :)
+    end function stdlib_linalg_${ri}$_expm_fun
+    #:endfor
+  end interface expm
+
+  !> Matrix exponential: subroutine interface
+  interface matrix_exp
+    !! version : experimental
+    !!
+    !! Computes the exponential of a matrix using a rational Pade approximation.
+    !! ([Specification](../page/specs/stdlib_linalg.html#matrix_exp))
+    !!
+    !! ### Description
+    !!
+    !! This interface provides methods for computing the exponential of a matrix
+    !! represented as a standard Fortran rank-2 array. Supported data types include
+    !! `real` and `complex`.
+    !!
+    !! By default, the order of the Pade approximation is set to 10. It can be changed
+    !! via the `order` argument that must be non-negative.
+    !!
+    !! If the input matrix is non-square or the order of the Pade approximation is
+    !! negative, the function returns an error state.
+    !!
+    !! ### Example
+    !!
+    !! ```fortran
+    !!  real(dp) :: A(3, 3), E(3, 3)
+    !!
+    !!  A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+    !!
+    !!  ! Default Pade approximation of the matrix exponential.
+    !!  call matrix_exp(A, E) ! Out-of-place
+    !!  ! call matrix_exp(A) for in-place computation.
+    !!
+    !!  ! Pade approximation with specified order.
+    !!  call matrix_exp(A, E, order=12)
+    !! ```
+    !!
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module subroutine stdlib_linalg_${ri}$_expm_inplace(A, order, err)
+        !> Input matrix A(n, n) / Output matrix E = exp(A)
+        ${rt}$, intent(inout) :: A(:, :)
+        !> [optional] Order of the Pade approximation (default `order=10`)
+        integer(ilp), optional, intent(in) :: order
+        !> [optional] Error handling.
+        type(linalg_state_type), optional, intent(out) :: err
+    end subroutine stdlib_linalg_${ri}$_expm_inplace
+
+    module subroutine stdlib_linalg_${ri}$_expm(A, E, order, err)
+        !> Input matrix A(n, n)
+        ${rt}$, intent(in) :: A(:, :)
+        !> Output matrix exponential E = exp(A)
+        ${rt}$, intent(out) :: E(:, :)
+        !> [optional] Order of the Pade approximation (default `order=10`)
+        integer(ilp), optional, intent(in) :: order
+        !> [optional] Error handling.
+        type(linalg_state_type), optional, intent(out) :: err
+    end subroutine stdlib_linalg_${ri}$_expm
+    #:endfor
+  end interface matrix_exp
 contains
 
 

src/stdlib_linalg_matrix_functions.fypp

@@ -0,0 +1,142 @@
+#:include "common.fypp"
+#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX, REAL_INIT))
+#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX, CMPLX_INIT))
+#:set RC_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_matrix_functions
+    use stdlib_constants
+    use stdlib_linalg_constants
+    use stdlib_linalg_blas, only: gemm
+    use stdlib_linalg_lapack, only: gesv, lacpy
+    use stdlib_linalg_lapack_aux, only: handle_gesv_info
+    use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
+    implicit none(type, external)
+
+    character(len=*), parameter :: this = "matrix_exponential"
+
+contains
+
+    #:for k,t,s, i in RC_KINDS_TYPES 
+    module function stdlib_linalg_${i}$_expm_fun(A, order) result(E)
+        !> Input matrix A(n, n).
+        ${t}$, intent(in) :: A(:, :)
+        !> [optional] Order of the Pade approximation.
+        integer(ilp), optional, intent(in) :: order
+        !> Exponential of the input matrix E = exp(A).
+        ${t}$, allocatable :: E(:, :)
+
+        E = A
+        call stdlib_linalg_${i}$_expm_inplace(E, order)
+    end function stdlib_linalg_${i}$_expm_fun
+
+    module subroutine stdlib_linalg_${i}$_expm(A, E, order, err)
+        !> Input matrix A(n, n).
+        ${t}$, intent(in) :: A(:, :)
+        !> Exponential of the input matrix E = exp(A).
+        ${t}$, intent(out) :: E(:, :)
+         !> [optional] Order of the Pade approximation.
+        integer(ilp), optional, intent(in) :: order
+        !> [optional] State return flag.
+        type(linalg_state_type), optional, intent(out) :: err
+
+        type(linalg_state_type) :: err0
+        integer(ilp) :: lda, n, lde, ne
+
+        ! Check E sizes
+        lda = size(A, 1, kind=ilp) ; n = size(A, 2, kind=ilp)
+        lde = size(E, 1, kind=ilp) ; ne = size(E, 2, kind=ilp)
+
+        if (lda<1 .or. n<1 .or. lda/=n .or. lde/=n .or. ne/=n) then     
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR, &
+                                     'invalid matrix sizes: A must be square (lda=', lda, ', n=', n, ')', &
+                                     ' E must be square (lde=', lde, ', ne=', ne, ')')
+        else
+            call lacpy("n", n, n, A, n, E, n) ! E = A
+            call stdlib_linalg_${i}$_expm_inplace(E, order, err0)
+        endif
+
+        ! Process output and return
+        call linalg_error_handling(err0,err)
+
+        return
+    end subroutine stdlib_linalg_${i}$_expm
+
+    module subroutine stdlib_linalg_${i}$_expm_inplace(A, order, err)
+        !> Input matrix A(n, n) / Output matrix exponential.
+        ${t}$, intent(inout) :: A(:, :)
+        !> [optional] Order of the Pade approximation.
+        integer(ilp), optional, intent(in) :: order
+        !> [optional] State return flag.
+        type(linalg_state_type), optional, intent(out) :: err
+
+        ! Internal variables.
+        ${t}$                   :: A2(size(A, 1), size(A, 2)), Q(size(A, 1), size(A, 2))
+        ${t}$                   :: X(size(A, 1), size(A, 2)), X_tmp(size(A, 1), size(A, 2))
+        real(${k}$)             :: a_norm, c
+        integer(ilp)            :: m, n, ee, k, s, order_, i, j
+        logical(lk)             :: p
+        type(linalg_state_type) :: err0
+
+        ! Deal with optional args.
+        order_ = 10 ; if (present(order)) order_ = order
+
+        ! Problem's dimension.
+        m = size(A, dim=1, kind=ilp) ; n = size(A, dim=2, kind=ilp)
+
+        if (m /= n) then
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'Invalid matrix size A=',[m, n])
+        else if (order_ < 0) then
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Order of Pade approximation &
+                                    needs to be positive, order=', order_)
+        else
+            ! Compute the L-infinity norm.
+            a_norm = mnorm(A, "inf")
+
+            ! Determine scaling factor for the matrix.
+            ee = int(log(a_norm) / log2_${k}$, kind=ilp) + 1
+            s  = max(0, ee+1)
+
+            ! Scale the input matrix & initialize polynomial.
+            A2 = A/2.0_${k}$**s
+            call lacpy("n", n, n, A2, n, X, n) ! X = A2
+
+            ! First step of the Pade approximation.
+            c = 0.5_${k}$
+            do concurrent(i=1:n, j=1:n)
+                A(i, j) = merge(1.0_${k}$ + c*A2(i, j), c*A2(i, j), i == j)
+                Q(i, j) = merge(1.0_${k}$ - c*A2(i, j), -c*A2(i, j), i == j)
+            enddo
+
+            ! Iteratively compute the Pade approximation.
+            p = .true.
+            do k = 2, order_
+                c = c * (order_ - k + 1) / (k * (2*order_ - k + 1))
+                call lacpy("n", n, n, X, n, X_tmp, n) ! X_tmp = X
+                call gemm("N", "N", n, n, n, one_${s}$, A2, n, X_tmp, n, zero_${s}$, X, n)
+                do concurrent(i=1:n, j=1:n)
+                    A(i, j) = A(i, j) + c*X(i, j)       ! E = E + c*X
+                    Q(i, j) = merge(Q(i, j) + c*X(i, j), Q(i, j) - c*X(i, j), p)
+                enddo
+                p = .not. p
+            enddo
+
+            block
+                integer(ilp) :: ipiv(n), info
+                call gesv(n, n, Q, n, ipiv, A, n, info) ! E = inv(Q) @ E
+                call handle_gesv_info(this, info, n, n, n, err0)
+            end block
+
+            ! Matrix squaring.
+            do k = 1, s
+                call lacpy("n", n, n, A, n, X, n) ! X = A
+                call gemm("N", "N", n, n, n, one_${s}$, X, n, X, n, zero_${s}$, A, n)
+            enddo
+        endif
+
+        call linalg_error_handling(err0, err)
+
+        return
+    end subroutine stdlib_linalg_${i}$_expm_inplace
+    #:endfor
+
+end submodule stdlib_linalg_matrix_functions