fix: exponential method ambiguity resolutions (#253)

* restrict types diagonal

* DefaultAlgorithm intercepts for Diagonal

* update to `default_algorithm`

* Algorithm selection

* Update src/implementations/exponential.jl

Co-authored-by: Lukas Devos <ldevos98@gmail.com>

* remove ambiguity

* Remove ambiguity

* alter ambiguity fix to be in line with decompositions

---------

Co-authored-by: Lukas Devos <ldevos98@gmail.com>

Author: sanderdemeyer

src/implementations/exponential.jl

@@ -8,14 +8,14 @@ copy_input(::typeof(exponential), A::Diagonal) = copy(A)
 copy_input(::typeof(exponential), (τ, A)::Tuple{Number, AbstractMatrix}) = (τ, copy!(similar(A, float(eltype(A))), A))
 copy_input(::typeof(exponential), (τ, A)::Tuple{Number, Diagonal}) = τ, copy(A)
 
-function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, alg::AbstractAlgorithm)
+function check_input(::typeof(exponential!), A::AbstractMatrix, expA, alg::AbstractAlgorithm)
     m = LinearAlgebra.checksquare(A)
     @check_size(expA, (m, m))
     @check_scalar(expA, A)
     return nothing
 end
 
-function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, ::DiagonalAlgorithm)
+function check_input(::typeof(exponential!), A::AbstractMatrix, expA, ::DiagonalAlgorithm)
     m = LinearAlgebra.checksquare(A)
     @assert isdiag(A)
     @assert expA isa Diagonal
@@ -24,14 +24,14 @@ function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMa
     return nothing
 end
 
-function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::AbstractAlgorithm)
+function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::AbstractAlgorithm)
     m = LinearAlgebra.checksquare(A)
     @check_size(expA, (m, m))
     @check_scalar(expA, A, (τ isa Real) ? identity : complex)
     return nothing
 end
 
-function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, ::DiagonalAlgorithm)
+function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, expA, ::DiagonalAlgorithm)
     m = LinearAlgebra.checksquare(A)
     @assert isdiag(A)
     @assert expA isa Diagonal
@@ -40,6 +40,13 @@ function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatr
     return nothing
 end
 
+# Algorithm selection
+# ---------------------------
+exponential!(A::AbstractMatrix, alg::DefaultAlgorithm) = exponential!(A, select_algorithm(exponential!, A, nothing; alg.kwargs...))
+exponential!(A::AbstractMatrix, out, alg::DefaultAlgorithm) = exponential!(A, out, select_algorithm(exponential!, A, nothing; alg.kwargs...))
+exponential!(τA::Tuple{Number, AbstractMatrix}, alg::DefaultAlgorithm) = exponential!(τA, select_algorithm(exponential!, τA, nothing; alg.kwargs...))
+exponential!(τA::Tuple{Number, AbstractMatrix}, out, alg::DefaultAlgorithm) = exponential!(τA, out, select_algorithm(exponential!, τA, nothing; alg.kwargs...))
+
 # Outputs
 # -------
 initialize_output(::typeof(exponential!), A::AbstractMatrix, ::AbstractAlgorithm) = A
@@ -48,22 +55,22 @@ initialize_output(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}
 
 # Implementation
 # --------------
-function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaLA)
+function exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaLA)
     check_input(exponential!, A, expA, alg)
     A = LinearAlgebra.exp!(A)
     A === expA || copy!(expA, A)
     return expA
 end
 
-exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEigh) = exponential!((one(eltype(A)), A), expA, alg)
-exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEig) = exponential!((one(eltype(A)), A), expA, alg)
+exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaEigh) = exponential!((one(eltype(A)), A), expA, alg)
+exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaEig) = exponential!((one(eltype(A)), A), expA, alg)
 
-function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::AbstractAlgorithm)
+function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::AbstractAlgorithm)
     expA .= A .* τ
     return exponential!(expA, expA, alg)
 end
 
-function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::MatrixFunctionViaEigh)
+function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::MatrixFunctionViaEigh)
     check_input(exponential!, (τ, A), expA, alg)
     D, V = eigh_full!(A, alg.eigh_alg)
     if eltype(A) <: Real
@@ -83,7 +90,7 @@ function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatr
     end
 end
 
-function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::MatrixFunctionViaEig)
+function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::MatrixFunctionViaEig)
     check_input(exponential!, (τ, A), expA, alg)
     D, V = eig_full!(A, alg.eig_alg)
     if eltype(A) <: Real && eltype(τ) <: Real
@@ -98,12 +105,12 @@ end
 
 # Diagonal logic
 # --------------
-function exponential!(A, expA, alg::DiagonalAlgorithm)
+function exponential!(A::AbstractMatrix, expA, alg::DiagonalAlgorithm)
     check_input(exponential!, A, expA, alg)
     return map_diagonal!(exp, expA, A)
 end
 
-function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::DiagonalAlgorithm)
+function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::DiagonalAlgorithm)
     check_input(exponential!, (τ, A), expA, alg)
     return map_diagonal!(x -> exp(x * τ), expA, A)
 end

src/interface/exponential.jl

@@ -35,5 +35,9 @@ function default_algorithm(::typeof(exponential!), ::Type{A}; kwargs...) where {
 end
 
 function default_algorithm(::typeof(exponential!), ::Tuple{A, B}; kwargs...) where {A, B}
-    return default_exponential_algorithm(B; kwargs...)
+    return default_algorithm(exponential!, B; kwargs...)
+end
+
+function default_algorithm(::typeof(exponential!), ::Type{Tuple{A, B}}; kwargs...) where {A, B}
+    return default_algorithm(exponential!, B; kwargs...)
 end