-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
8 changed files
with
135 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,17 @@ | ||
# Expm | ||
|
||
This section of the manual describes the `Expm` module. | ||
|
||
```@contents | ||
Pages = ["Expm.md"] | ||
Depth = 3 | ||
``` | ||
|
||
```@meta | ||
CurrentModule = ReachabilityBase.Expm | ||
``` | ||
|
||
```@docs | ||
Expm | ||
expm | ||
``` |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,82 @@ | ||
""" | ||
Expm | ||
This module provides functionality for computing the matrix exponential. | ||
""" | ||
module Expm | ||
|
||
using Requires | ||
using LinearAlgebra: checksquare | ||
using SparseArrays: AbstractSparseMatrixCSC | ||
|
||
export expm | ||
|
||
# matrix dimension from which on the FastExpm algorithm is used | ||
const FASTEXPM_THRESHOLD = 17 | ||
|
||
""" | ||
expm(M::AbstractMatrix) | ||
Compute the matrix exponential. Conditionally (see below) use the | ||
[`FastExpm`](https://github.com/fmentink/FastExpm.jl) implementation. | ||
### Input | ||
- `M` -- matrix | ||
### Output | ||
A matrix corresponding to ``\\exp(M)``. | ||
### Algorithm | ||
The [`FastExpm`](https://github.com/fmentink/FastExpm.jl) implementation is | ||
used under the following conditions: | ||
- The `FastExpm` package is loaded. | ||
- The matrix `M` is sparse. | ||
- The dimension of `M` is higher than $FASTEXPM_THRESHOLD. | ||
### Notes | ||
If not all of the above conditions are satisfied, this method uses the | ||
`LinearAlgebra` implementation, for which it converts to a dense matrix. | ||
If the `FastExpm` implementation is used, the result is a complex matrix. If | ||
`M` is a real matrix, the result is a real matrix too (i.e., the imaginary | ||
part is discarded). | ||
""" | ||
function expm(M::AbstractMatrix) | ||
return exp(Matrix(M)) # convert to dense matrix | ||
end | ||
|
||
# skip conversion to dense matrix if already dense | ||
function expm(M::Matrix) | ||
return exp(M) | ||
end | ||
|
||
function load_FastExpm() | ||
return quote | ||
using .FastExpm: fastExpm | ||
function expm(M::AbstractSparseMatrixCSC) | ||
if checksquare(M) >= FASTEXPM_THRESHOLD | ||
return fastExpm(M) | ||
else | ||
return exp(Matrix(M)) # convert to dense matrix | ||
end | ||
end | ||
function expm(M::AbstractSparseMatrixCSC{<:Real}) | ||
if checksquare(M) >= FASTEXPM_THRESHOLD | ||
return real.(fastExpm(M)) # drop imaginary part | ||
else | ||
return exp(Matrix(M)) # convert to dense matrix | ||
end | ||
end | ||
end | ||
end | ||
|
||
function __init__() | ||
@require FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" eval(load_FastExpm()) | ||
end | ||
|
||
end # module |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,27 @@ | ||
using ReachabilityBase.Expm | ||
using SparseArrays | ||
|
||
M1 = [1.0 0; 0 3] # small matrix | ||
M2 = [1.0+1im 0; 0 3] # complex | ||
M3 = zeros(20, 20) | ||
for i in 1:20 # large matrix | ||
M3[i, i] = 1 / i | ||
end | ||
M4 = fill(0.0+0im, 20, 20) # complex | ||
for i in 1:20 | ||
M4[i, i] = 1 / i | ||
end | ||
|
||
for M in (M1, M2, M3, M4) | ||
@test expm(M) == expm(sparse(M)) == exp(M) | ||
end | ||
|
||
# switch to FastExpm | ||
import FastExpm | ||
|
||
for M in (M1, M2) | ||
@test expm(M) == expm(sparse(M)) == exp(M) | ||
end | ||
for M in (M3, M4) | ||
@test expm(M) ≈ expm(sparse(M)) ≈ exp(M) | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,9 +1,11 @@ | ||
[deps] | ||
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" | ||
FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" | ||
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" | ||
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" | ||
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" | ||
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" | ||
|
||
[compat] | ||
Documenter = "0.27" | ||
FastExpm = "1.1" |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters