-
Notifications
You must be signed in to change notification settings - Fork 191
Exponential of a matrix #968
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: master
Are you sure you want to change the base?
Changes from 4 commits
b0a74b1
fa36f33
d8b1857
4089d18
c6857bc
4310db5
59ffb20
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -52,3 +52,4 @@ ADD_EXAMPLE(qr) | |
ADD_EXAMPLE(qr_space) | ||
ADD_EXAMPLE(cholesky) | ||
ADD_EXAMPLE(chol) | ||
ADD_EXAMPLE(expm) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,7 @@ | ||
program example_expm | ||
use stdlib_linalg, only: expm | ||
implicit none | ||
real :: A(3, 3), E(3, 3) | ||
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) | ||
E = expm(A) | ||
end program example_expm |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,130 @@ | ||
#: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_lapack, only: gesv | ||
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & | ||
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR | ||
implicit none | ||
|
||
#:for rk, rt, ri in (REAL_KINDS_TYPES) | ||
${rt}$, parameter :: zero_${ri}$ = 0._${rk}$ | ||
${rt}$, parameter :: one_${ri}$ = 1._${rk}$ | ||
#:endfor | ||
#:for rk, rt, ri in (CMPLX_KINDS_TYPES) | ||
${rt}$, parameter :: zero_${ri}$ = (0._${rk}$, 0._${rk}$) | ||
${rt}$, parameter :: one_${ri}$ = (1._${rk}$, 0._${rk}$) | ||
#:endfor | ||
|
||
contains | ||
|
||
#:for rk,rt,ri in RC_KINDS_TYPES | ||
module function stdlib_expm_${ri}$(A, order, err) result(E) | ||
!> Input matrix A(n, n). | ||
${rt}$, intent(in) :: 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 | ||
!> Exponential of the input matrix E = exp(A). | ||
${rt}$, allocatable :: E(:, :) | ||
|
||
! Internal variables. | ||
${rt}$, allocatable :: A2(:, :), Q(:, :), X(:, :) | ||
real(${rk}$) :: a_norm, c | ||
integer(ilp) :: m, n, ee, k, s, order_, i, j | ||
logical(lk) :: p | ||
character(len=*), parameter :: this = "expm" | ||
type(linalg_state_type) :: err0 | ||
|
||
! Deal with optional args. | ||
order_ = 10 ; if (present(order)) order_ = order | ||
|
||
! Problem's dimension. | ||
m = size(A, 1) ; n = size(A, 2) | ||
|
||
if (m /= n) then | ||
err = linalg_state_type(this,LINALG_VALUE_ERROR,'Invalid matrix size A=',[m, n]) | ||
call linalg_error_handling(err0, err) | ||
else if (order_ < 0) then | ||
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Order of Pade approximation & | ||
needs to be positive, order=', order_) | ||
call linalg_error_handling(err0, err) | ||
endif | ||
|
||
! 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 | ||
|
||
! First step of the Pade approximation. | ||
c = 0.5_${rk}$ | ||
allocate (E, source=A2) ; allocate (Q, source=A2) | ||
do concurrent(i=1:n, j=1:n) | ||
E(i, j) = c*E(i, j) ; if (i == j) E(i, j) = 1.0_${rk}$ + E(i, j) ! E = I + c*A2 | ||
Q(i, j) = -c*Q(i, j) ; if (i == j) Q(i, j) = 1.0_${rk}$ + Q(i, j) ! Q = I - c*A2 | ||
enddo | ||
|
||
! Iteratively compute the Pade approximation. | ||
p = .true. | ||
do k = 2, order_ | ||
c = c * (order_ - k + 1) / (k * (2*order_ - k + 1)) | ||
X = matmul(A2, X) | ||
do concurrent(i=1:n, j=1:n) | ||
E(i, j) = E(i, j) + c*X(i, j) ! E = E + c*X | ||
enddo | ||
if (p) then | ||
do concurrent(i=1:n, j=1:n) | ||
Q(i, j) = Q(i, j) + c*X(i, j) ! Q = Q + c*X | ||
enddo | ||
else | ||
do concurrent(i=1:n, j=1:n) | ||
Q(i, j) = Q(i, j) - c*X(i, j) ! Q = Q - c*X | ||
enddo | ||
endif | ||
p = .not. p | ||
enddo | ||
|
||
block | ||
integer(ilp) :: ipiv(n), info | ||
call gesv(n, n, Q, n, ipiv, E, n, info) ! E = inv(Q) @ E | ||
call handle_gesv_info(info, n, n, n, err0) | ||
call linalg_error_handling(err0, err) | ||
end block | ||
|
||
! This loop should eventually be replaced by a fast matrix_power function. | ||
do k = 1, s | ||
E = matmul(E, E) | ||
loiseaujc marked this conversation as resolved.
Show resolved
Hide resolved
|
||
enddo | ||
return | ||
contains | ||
elemental subroutine handle_gesv_info(info,lda,n,nrhs,err) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should this procedure not provided by There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think it should. I've sent a PM on the discourse to discuss this issue actually. I was suggesting to create a new There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, I agree that code duplication should be avoided. Looking at the current lapack module structure, there are so many modules and submodules already. So I think that we should try to not create a new one only for the error message handling? Perhaps we should use |
||
integer(ilp), intent(in) :: info,lda,n,nrhs | ||
type(linalg_state_type), intent(out) :: err | ||
! Process output | ||
select case (info) | ||
case (0) | ||
! Success | ||
case (-1) | ||
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n) | ||
case (-2) | ||
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs) | ||
case (-4) | ||
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n]) | ||
case (-7) | ||
err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n]) | ||
case (1:) | ||
err = linalg_state_type(this,LINALG_ERROR,'singular matrix') | ||
case default | ||
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') | ||
end select | ||
end subroutine handle_gesv_info | ||
end function stdlib_expm_${ri}$ | ||
#:endfor | ||
|
||
end submodule stdlib_linalg_matrix_functions |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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) | ||
loiseaujc marked this conversation as resolved.
Show resolved
Hide resolved
|
||
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 |
Uh oh!
There was an error while loading. Please reload this page.