Skip to content

Commit 56474e1

Browse files
loiseaujcloiseaujc
authored
Merge branch 'matrix_exponential' into master
2 parents 089d316 + 59ffb20 commit 56474e1

File tree

8 files changed

+361
-3
lines changed

8 files changed

+361
-3
lines changed

‎doc/specs/stdlib_linalg.md

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1884,3 +1884,38 @@ 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 compute its matrix exponential \(E = \exp(A)\) using a Pade approximation.
1896+
1897+
### Syntax
1898+
1899+
`E = ` [[stdlib_linalg(module):expm(interface)]] `(a [, order, err])`
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+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
1908+
1909+
### Return value
1910+
1911+
The returned array `E` contains the Pade approximation of \(\exp(A)\).
1912+
1913+
If `A` is non-square or `order` is negative, it raise a `LINALG_VALUE_ERROR`.
1914+
If `err` is not present, exceptions trigger an `error stop`.
1915+
1916+
### Example
1917+
1918+
```fortran
1919+
{!example/linalg/example_expm.f90!}
1920+
```
1921+

‎example/linalg/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,3 +52,4 @@ ADD_EXAMPLE(qr)
5252
ADD_EXAMPLE(qr_space)
5353
ADD_EXAMPLE(cholesky)
5454
ADD_EXAMPLE(chol)
55+
ADD_EXAMPLE(expm)

‎example/linalg/example_expm.f90

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
program example_expm
2+
use stdlib_linalg, only: expm
3+
implicit none
4+
real :: A(3, 3), E(3, 3)
5+
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
6+
E = expm(A)
7+
end program example_expm

‎src/CMakeLists.txt

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,16 +32,17 @@ set(fppFiles
3232
stdlib_linalg_kronecker.fypp
3333
stdlib_linalg_cross_product.fypp
3434
stdlib_linalg_eigenvalues.fypp
35-
stdlib_linalg_solve.fypp
35+
stdlib_linalg_solve.fypp
3636
stdlib_linalg_determinant.fypp
3737
stdlib_linalg_qr.fypp
3838
stdlib_linalg_inverse.fypp
3939
stdlib_linalg_pinv.fypp
4040
stdlib_linalg_norms.fypp
4141
stdlib_linalg_state.fypp
42-
stdlib_linalg_svd.fypp
42+
stdlib_linalg_svd.fypp
4343
stdlib_linalg_cholesky.fypp
4444
stdlib_linalg_schur.fypp
45+
stdlib_linalg_matrix_functions.fypp
4546
stdlib_optval.fypp
4647
stdlib_selection.fypp
4748
stdlib_sorting.fypp

‎src/stdlib_linalg.fypp

Lines changed: 48 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
3132
public :: eye
3233
public :: inv
3334
public :: invert
@@ -1678,6 +1679,53 @@ 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 which 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_expm_${ri}$(A, order, err) result(E)
1717+
!> Input matrix a(n, n).
1718+
${rt}$, intent(in) :: A(:, :)
1719+
!> [optional] Order of the Pade approximation (default `order=10`)
1720+
integer(ilp), optional, intent(in) :: order
1721+
!> [optional] State return flag. On error, if not requested, the code will stop.
1722+
type(linalg_state_type), optional, intent(out) :: err
1723+
!> Exponential of the input matrix E = exp(A).
1724+
${rt}$, allocatable :: E(:, :)
1725+
end function stdlib_expm_${ri}$
1726+
#:endfor
1727+
end interface expm
1728+
16811729
contains
16821730

16831731

‎src/stdlib_linalg_matrix_functions.fypp

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

‎test/linalg/CMakeLists.txt

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ set(
1515
"test_linalg_matrix_property_checks.fypp"
1616
"test_linalg_sparse.fypp"
1717
"test_linalg_cholesky.fypp"
18+
"test_linalg_expm.fypp"
1819
)
1920

2021
# Preprocessed files to contain preprocessor directives -> .F90
@@ -41,4 +42,5 @@ ADDTEST(linalg_qr)
4142
ADDTEST(linalg_schur)
4243
ADDTEST(linalg_svd)
4344
ADDTEST(linalg_sparse)
44-
ADDTESTPP(blas_lapack)
45+
ADDTEST(linalg_expm)
46+
ADDTESTPP(blas_lapack)

0 commit comments

Comments
 (0)