Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Sign up
Appearance settings

Commit 7ee5aa8

Browse files
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
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

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)

AltStyle によって変換されたページ (->オリジナル) /