Taylor series native exponential implementation

benchmark/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"

benchmark/README.md

@@ -0,0 +1,47 @@
+# Matrix exponential benchmarks
+
+`bench_exponential.jl` times and accuracy-tests every applicable `MatrixAlgebraKit.exponential`
+algorithm (`MatrixFunctionViaTaylor`, `MatrixFunctionViaLA`, `MatrixFunctionViaEig`,
+`MatrixFunctionViaEigh`) for `Float64`, `ComplexF64`, `Double64`, `Complex{Double64}`,
+`BigFloat` and `Complex{BigFloat}`.
+The `Double64` (double-double, ~106-bit) types sit between the hardware BLAS floats and
+arbitrary-precision `BigFloat`, and exercise the `Eig`/`Eigh` paths through the generic
+GenericSchur / GenericLinearAlgebra drivers (the `LA` LAPACK path applies to BLAS types only).
+
+Accuracy inputs are built from a *known* spectrum as `A = V · Diag(λ) · V⁻¹`, so the exact
+`exp(A) = V · Diag(exp λ) · V⁻¹` is available analytically (computed at high precision).
+Eigenvalue families (`small`, `wide`, `imaginary`, `stiff`, `illconditioned`, `hermitian`)
+and eigenvector families (`unitary`, `general`, `illconditioned`) are chosen to stress or
+favor particular algorithms.
+
+## Setup (once)
+
+From the repository root:
+
+```sh
+julia --project=benchmark -e 'using Pkg; Pkg.develop(path=pwd()); Pkg.instantiate()'
+```
+
+This `dev`s the in-tree `MatrixAlgebraKit` into the benchmark's own environment.
+
+## Run
+
+```sh
+julia --project=benchmark benchmark/bench_exponential.jl            # full run
+julia --project=benchmark benchmark/bench_exponential.jl --quick    # fast smoke run
+julia --project=benchmark benchmark/bench_exponential.jl --csv out.csv
+```
+
+Results are printed as per-type accuracy and timing tables and written to two CSVs derived
+from the `--csv` path (default base `benchmark/results.csv`): `*-accuracy.csv` and
+`*-timing.csv`. They are separate files because the accuracy and timing passes report
+disjoint quantities — a single combined table would leave the timing columns empty on the
+accuracy rows and vice-versa.
+
+## Accuracy columns
+
+- `analytic_err` — relative error vs the analytic `V·Diag(exp λ)·V⁻¹` reference (primary).
+- `recomp_err` — relative error vs an independent high-precision `exp` of the exact input.
+- `expexp_err` — `‖exp(A)·exp(−A) − I‖₁` (reference-free consistency check).
+- `det_err` — relative error of `det(exp A)` against `exp(tr A)`.
+- `condV` — designed condition number of the eigenvector basis (context for non-normal cases).