Add Expm module

.github/workflows/ci_PR.yml

@@ -22,7 +22,8 @@ jobs:
         arch:
           - x64
         include:
-          - version: '1.2'  # test on oldest supported version
+          #- version: '1.2'  # test on oldest supported version
+          - version: '1.5'  # test on oldest version supporting all test dependencies
             arch: x64
             os: ubuntu-latest
     env:

docs/make.jl

@@ -20,7 +20,8 @@ makedocs(; sitename="ReachabilityBase.jl",
                                  "Distribution" => "lib/Distribution.md",
                                  "Subtypes" => "lib/Subtypes.md",
                                  "Arrays" => "lib/Arrays.md",
-                                 "Timing" => "lib/Timing.md"],
+                                 "Timing" => "lib/Timing.md",
+                                 "Expm" => "lib/Expm.md"],
                 "About" => "about.md"],
          doctest=false,
          strict=true)

docs/src/lib/Expm.md

@@ -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
+```

src/Expm/Expm.jl

@@ -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

src/ReachabilityBase.jl

@@ -9,5 +9,6 @@ include("Distribution/Distribution.jl")
 include("Subtypes/Subtypes.jl")
 include("Arrays/Arrays.jl")
 include("Timing/Timing.jl")
+include("Expm/Expm.jl")
 
 end  # module

test/Expm/expm.jl

@@ -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

test/Project.toml

@@ -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"

test/runtests.jl

@@ -28,6 +28,9 @@ end
 @time @testset "Arrays.vector_operations" begin
     include("Arrays/vector_operations.jl")
 end
+@time @testset "Expm.expm" begin
+    include("Expm/expm.jl")
+end
 
 if VERSION > v"1.6"
     using Documenter