Skip to content

Commit 7ee5aa8

Browse files
jvdp1jvdp1
authored
Addition of matrix exponential
Matrix exponential
2 parents e5d296d + 305067c commit 7ee5aa8

File tree

10 files changed

+516
-0
lines changed

10 files changed

+516
-0
lines changed

‎doc/specs/stdlib_linalg.md‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1884,3 +1884,72 @@ If `err` is not present, exceptions trigger an `error stop`.
18841884
{!example/linalg/example_mnorm.f90!}
18851885
```
18861886

1887+
## `expm` - Computes the matrix exponential {#expm}
1888+
1889+
### Status
1890+
1891+
Experimental
1892+
1893+
### Description
1894+
1895+
Given a matrix \(A\), this function computes its matrix exponential \(E = \exp(A)\) using a Pade approximation.
1896+
1897+
### Syntax
1898+
1899+
`E = ` [[stdlib_linalg(module):expm(interface)]] `(a [, order])`
1900+
1901+
### Arguments
1902+
1903+
`a`: Shall be a rank-2 `real` or `complex` array containing the data. It is an `intent(in)` argument.
1904+
1905+
`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.
1906+
1907+
### Return value
1908+
1909+
The returned array `E` contains the Pade approximation of \(\exp(A)\).
1910+
1911+
If `A` is non-square or `order` is negative, it raises a `LINALG_VALUE_ERROR`.
1912+
1913+
### Example
1914+
1915+
```fortran
1916+
{!example/linalg/example_expm.f90!}
1917+
```
1918+
1919+
## `matrix_exp` - Computes the matrix exponential {#matrix_exp}
1920+
1921+
### Status
1922+
1923+
Experimental
1924+
1925+
### Description
1926+
1927+
Given a matrix \(A\), this function computes its matrix exponential \(E = \exp(A)\) using a Pade approximation.
1928+
1929+
### Syntax
1930+
1931+
`call ` [[stdlib_linalg(module):matrix_exp(interface)]] `(a [, e, order, err])`
1932+
1933+
### Arguments
1934+
1935+
`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.
1936+
1937+
`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`.
1938+
1939+
`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.
1940+
1941+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
1942+
1943+
### Return value
1944+
1945+
The returned array `A` (in-place) or `E` (out-of-place) contains the Pade approximation of \(\exp(A)\).
1946+
1947+
If `A` is non-square or `order` is negative, it raises a `LINALG_VALUE_ERROR`.
1948+
If `err` is not present, exceptions trigger an `error stop`.
1949+
1950+
### Example
1951+
1952+
```fortran
1953+
{!example/linalg/example_matrix_exp.f90!}
1954+
```
1955+

‎example/linalg/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,3 +57,5 @@ ADD_EXAMPLE(qr)
5757
ADD_EXAMPLE(qr_space)
5858
ADD_EXAMPLE(cholesky)
5959
ADD_EXAMPLE(chol)
60+
ADD_EXAMPLE(expm)
61+
ADD_EXAMPLE(matrix_exp)

‎example/linalg/example_expm.f90‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
program example_expm
2+
use stdlib_linalg, only: expm
3+
implicit none
4+
real :: A(3, 3), E(3, 3)
5+
integer :: i
6+
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
7+
E = expm(A)
8+
9+
print *, "Matrix A :"
10+
do i = 1, 3
11+
print *, A(i, :)
12+
end do
13+
14+
print *, "Matrix exponential E = exp(A):"
15+
do i = 1, 3
16+
print *, E(i, :)
17+
end do
18+
end program example_expm

‎example/linalg/example_matrix_exp.f90‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
program example_expm
2+
use stdlib_linalg, only: matrix_exp
3+
implicit none
4+
real :: A(3, 3), E(3, 3)
5+
integer :: i
6+
7+
print *, "Matrix A :"
8+
do i = 1, 3
9+
print *, A(i, :)
10+
end do
11+
12+
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
13+
call matrix_exp(A) ! In-place computation.
14+
! For out-of-place, use call matrix_exp(A, E).
15+
16+
print *, "Matrix exponential E = exp(A):"
17+
do i = 1, 3
18+
print *, E(i, :)
19+
end do
20+
end program example_expm

‎src/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@ set(fppFiles
4646
stdlib_linalg_svd.fypp
4747
stdlib_linalg_cholesky.fypp
4848
stdlib_linalg_schur.fypp
49+
stdlib_linalg_matrix_functions.fypp
4950
stdlib_optval.fypp
5051
stdlib_selection.fypp
5152
stdlib_sorting.fypp

‎src/stdlib_constants.fypp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,7 @@ module stdlib_constants
7171
#:for k, t, s in R_KINDS_TYPES
7272
${t}$, parameter, public :: zero_${s}$ = 0._${k}$
7373
${t}$, parameter, public :: one_${s}$ = 1._${k}$
74+
${t}$, parameter, public :: log2_${s}$ = log(2.0_${k}$)
7475
#:endfor
7576
#:for k, t, s in C_KINDS_TYPES
7677
${t}$, parameter, public :: zero_${s}$ = (0._${k}$,0._${k}$)

‎src/stdlib_linalg.fypp‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ module stdlib_linalg
2828
public :: eigh
2929
public :: eigvals
3030
public :: eigvalsh
31+
public :: expm, matrix_exp
3132
public :: eye
3233
public :: inv
3334
public :: invert
@@ -1678,6 +1679,107 @@ module stdlib_linalg
16781679
#:endfor
16791680
end interface mnorm
16801681

1682+
!> Matrix exponential: function interface
1683+
interface expm
1684+
!! version : experimental
1685+
!!
1686+
!! Computes the exponential of a matrix using a rational Pade approximation.
1687+
!! ([Specification](../page/specs/stdlib_linalg.html#expm))
1688+
!!
1689+
!! ### Description
1690+
!!
1691+
!! This interface provides methods for computing the exponential of a matrix
1692+
!! represented as a standard Fortran rank-2 array. Supported data types include
1693+
!! `real` and `complex`.
1694+
!!
1695+
!! By default, the order of the Pade approximation is set to 10. It can be changed
1696+
!! via the `order` argument that must be non-negative.
1697+
!!
1698+
!! If the input matrix is non-square or the order of the Pade approximation is
1699+
!! negative, the function returns an error state.
1700+
!!
1701+
!! ### Example
1702+
!!
1703+
!! ```fortran
1704+
!! real(dp) :: A(3, 3), E(3, 3)
1705+
!!
1706+
!! A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
1707+
!!
1708+
!! ! Default Pade approximation of the matrix exponential.
1709+
!! E = expm(A)
1710+
!!
1711+
!! ! Pade approximation with specified order.
1712+
!! E = expm(A, order=12)
1713+
!! ```
1714+
!!
1715+
#:for rk,rt,ri in RC_KINDS_TYPES
1716+
module function stdlib_linalg_${ri}$_expm_fun(A, order) result(E)
1717+
!> Input matrix a(:, :).
1718+
${rt}$, intent(in) :: A(:, :)
1719+
!> [optional] Order of the Pade approximation (default `order=10`)
1720+
integer(ilp), optional, intent(in) :: order
1721+
!> Exponential of the input matrix E = exp(A).
1722+
${rt}$, allocatable :: E(:, :)
1723+
end function stdlib_linalg_${ri}$_expm_fun
1724+
#:endfor
1725+
end interface expm
1726+
1727+
!> Matrix exponential: subroutine interface
1728+
interface matrix_exp
1729+
!! version : experimental
1730+
!!
1731+
!! Computes the exponential of a matrix using a rational Pade approximation.
1732+
!! ([Specification](../page/specs/stdlib_linalg.html#matrix_exp))
1733+
!!
1734+
!! ### Description
1735+
!!
1736+
!! This interface provides methods for computing the exponential of a matrix
1737+
!! represented as a standard Fortran rank-2 array. Supported data types include
1738+
!! `real` and `complex`.
1739+
!!
1740+
!! By default, the order of the Pade approximation is set to 10. It can be changed
1741+
!! via the `order` argument that must be non-negative.
1742+
!!
1743+
!! If the input matrix is non-square or the order of the Pade approximation is
1744+
!! negative, the function returns an error state.
1745+
!!
1746+
!! ### Example
1747+
!!
1748+
!! ```fortran
1749+
!! real(dp) :: A(3, 3), E(3, 3)
1750+
!!
1751+
!! A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
1752+
!!
1753+
!! ! Default Pade approximation of the matrix exponential.
1754+
!! call matrix_exp(A, E) ! Out-of-place
1755+
!! ! call matrix_exp(A) for in-place computation.
1756+
!!
1757+
!! ! Pade approximation with specified order.
1758+
!! call matrix_exp(A, E, order=12)
1759+
!! ```
1760+
!!
1761+
#:for rk,rt,ri in RC_KINDS_TYPES
1762+
module subroutine stdlib_linalg_${ri}$_expm_inplace(A, order, err)
1763+
!> Input matrix A(n, n) / Output matrix E = exp(A)
1764+
${rt}$, intent(inout) :: A(:, :)
1765+
!> [optional] Order of the Pade approximation (default `order=10`)
1766+
integer(ilp), optional, intent(in) :: order
1767+
!> [optional] Error handling.
1768+
type(linalg_state_type), optional, intent(out) :: err
1769+
end subroutine stdlib_linalg_${ri}$_expm_inplace
1770+
1771+
module subroutine stdlib_linalg_${ri}$_expm(A, E, order, err)
1772+
!> Input matrix A(n, n)
1773+
${rt}$, intent(in) :: A(:, :)
1774+
!> Output matrix exponential E = exp(A)
1775+
${rt}$, intent(out) :: E(:, :)
1776+
!> [optional] Order of the Pade approximation (default `order=10`)
1777+
integer(ilp), optional, intent(in) :: order
1778+
!> [optional] Error handling.
1779+
type(linalg_state_type), optional, intent(out) :: err
1780+
end subroutine stdlib_linalg_${ri}$_expm
1781+
#:endfor
1782+
end interface matrix_exp
16811783
contains
16821784

16831785

‎src/stdlib_linalg_matrix_functions.fypp‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX, REAL_INIT))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX, CMPLX_INIT))
4+
#:set RC_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES
5+
submodule (stdlib_linalg) stdlib_linalg_matrix_functions
6+
use stdlib_constants
7+
use stdlib_linalg_constants
8+
use stdlib_linalg_blas, only: gemm
9+
use stdlib_linalg_lapack, only: gesv, lacpy
10+
use stdlib_linalg_lapack_aux, only: handle_gesv_info
11+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
12+
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
13+
implicit none(type, external)
14+
15+
character(len=*), parameter :: this = "matrix_exponential"
16+
17+
contains
18+
19+
#:for k,t,s, i in RC_KINDS_TYPES
20+
module function stdlib_linalg_${i}$_expm_fun(A, order) result(E)
21+
!> Input matrix A(n, n).
22+
${t}$, intent(in) :: A(:, :)
23+
!> [optional] Order of the Pade approximation.
24+
integer(ilp), optional, intent(in) :: order
25+
!> Exponential of the input matrix E = exp(A).
26+
${t}$, allocatable :: E(:, :)
27+
28+
E = A
29+
call stdlib_linalg_${i}$_expm_inplace(E, order)
30+
end function stdlib_linalg_${i}$_expm_fun
31+
32+
module subroutine stdlib_linalg_${i}$_expm(A, E, order, err)
33+
!> Input matrix A(n, n).
34+
${t}$, intent(in) :: A(:, :)
35+
!> Exponential of the input matrix E = exp(A).
36+
${t}$, intent(out) :: E(:, :)
37+
!> [optional] Order of the Pade approximation.
38+
integer(ilp), optional, intent(in) :: order
39+
!> [optional] State return flag.
40+
type(linalg_state_type), optional, intent(out) :: err
41+
42+
type(linalg_state_type) :: err0
43+
integer(ilp) :: lda, n, lde, ne
44+
45+
! Check E sizes
46+
lda = size(A, 1, kind=ilp) ; n = size(A, 2, kind=ilp)
47+
lde = size(E, 1, kind=ilp) ; ne = size(E, 2, kind=ilp)
48+
49+
if (lda<1 .or. n<1 .or. lda/=n .or. lde/=n .or. ne/=n) then
50+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR, &
51+
'invalid matrix sizes: A must be square (lda=', lda, ', n=', n, ')', &
52+
' E must be square (lde=', lde, ', ne=', ne, ')')
53+
else
54+
call lacpy("n", n, n, A, n, E, n) ! E = A
55+
call stdlib_linalg_${i}$_expm_inplace(E, order, err0)
56+
endif
57+
58+
! Process output and return
59+
call linalg_error_handling(err0,err)
60+
61+
return
62+
end subroutine stdlib_linalg_${i}$_expm
63+
64+
module subroutine stdlib_linalg_${i}$_expm_inplace(A, order, err)
65+
!> Input matrix A(n, n) / Output matrix exponential.
66+
${t}$, intent(inout) :: A(:, :)
67+
!> [optional] Order of the Pade approximation.
68+
integer(ilp), optional, intent(in) :: order
69+
!> [optional] State return flag.
70+
type(linalg_state_type), optional, intent(out) :: err
71+
72+
! Internal variables.
73+
${t}$ :: A2(size(A, 1), size(A, 2)), Q(size(A, 1), size(A, 2))
74+
${t}$ :: X(size(A, 1), size(A, 2)), X_tmp(size(A, 1), size(A, 2))
75+
real(${k}$) :: a_norm, c
76+
integer(ilp) :: m, n, ee, k, s, order_, i, j
77+
logical(lk) :: p
78+
type(linalg_state_type) :: err0
79+
80+
! Deal with optional args.
81+
order_ = 10 ; if (present(order)) order_ = order
82+
83+
! Problem's dimension.
84+
m = size(A, dim=1, kind=ilp) ; n = size(A, dim=2, kind=ilp)
85+
86+
if (m /= n) then
87+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'Invalid matrix size A=',[m, n])
88+
else if (order_ < 0) then
89+
err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Order of Pade approximation &
90+
needs to be positive, order=', order_)
91+
else
92+
! Compute the L-infinity norm.
93+
a_norm = mnorm(A, "inf")
94+
95+
! Determine scaling factor for the matrix.
96+
ee = int(log(a_norm) / log2_${k}$, kind=ilp) + 1
97+
s = max(0, ee+1)
98+
99+
! Scale the input matrix & initialize polynomial.
100+
A2 = A/2.0_${k}$**s
101+
call lacpy("n", n, n, A2, n, X, n) ! X = A2
102+
103+
! First step of the Pade approximation.
104+
c = 0.5_${k}$
105+
do concurrent(i=1:n, j=1:n)
106+
A(i, j) = merge(1.0_${k}$ + c*A2(i, j), c*A2(i, j), i == j)
107+
Q(i, j) = merge(1.0_${k}$ - c*A2(i, j), -c*A2(i, j), i == j)
108+
enddo
109+
110+
! Iteratively compute the Pade approximation.
111+
p = .true.
112+
do k = 2, order_
113+
c = c * (order_ - k + 1) / (k * (2*order_ - k + 1))
114+
call lacpy("n", n, n, X, n, X_tmp, n) ! X_tmp = X
115+
call gemm("N", "N", n, n, n, one_${s}$, A2, n, X_tmp, n, zero_${s}$, X, n)
116+
do concurrent(i=1:n, j=1:n)
117+
A(i, j) = A(i, j) + c*X(i, j) ! E = E + c*X
118+
Q(i, j) = merge(Q(i, j) + c*X(i, j), Q(i, j) - c*X(i, j), p)
119+
enddo
120+
p = .not. p
121+
enddo
122+
123+
block
124+
integer(ilp) :: ipiv(n), info
125+
call gesv(n, n, Q, n, ipiv, A, n, info) ! E = inv(Q) @ E
126+
call handle_gesv_info(this, info, n, n, n, err0)
127+
end block
128+
129+
! Matrix squaring.
130+
do k = 1, s
131+
call lacpy("n", n, n, A, n, X, n) ! X = A
132+
call gemm("N", "N", n, n, n, one_${s}$, X, n, X, n, zero_${s}$, A, n)
133+
enddo
134+
endif
135+
136+
call linalg_error_handling(err0, err)
137+
138+
return
139+
end subroutine stdlib_linalg_${i}$_expm_inplace
140+
#:endfor
141+
142+
end submodule stdlib_linalg_matrix_functions

0 commit comments

Comments
 (0)