|
| 1 | +# Inputs |
| 2 | +# ------ |
| 3 | +function copy_input(::typeof(exponential), A::AbstractMatrix) |
| 4 | + return copy!(similar(A, float(eltype(A))), A) |
| 5 | +end |
| 6 | + |
| 7 | +copy_input(::typeof(exponential), A::Diagonal) = copy(A) |
| 8 | +copy_input(::typeof(exponential), (τ, A)::Tuple{Number, AbstractMatrix}) = (τ, copy!(similar(A, float(eltype(A))), A)) |
| 9 | +copy_input(::typeof(exponential), (τ, A)::Tuple{Number, Diagonal}) = τ, copy(A) |
| 10 | + |
| 11 | +function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, alg::AbstractAlgorithm) |
| 12 | + m = LinearAlgebra.checksquare(A) |
| 13 | + @check_size(expA, (m, m)) |
| 14 | + @check_scalar(expA, A) |
| 15 | + return nothing |
| 16 | +end |
| 17 | + |
| 18 | +function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, ::DiagonalAlgorithm) |
| 19 | + m = LinearAlgebra.checksquare(A) |
| 20 | + @assert isdiag(A) |
| 21 | + @assert expA isa Diagonal |
| 22 | + @check_size(expA, (m, m)) |
| 23 | + @check_scalar(expA, A) |
| 24 | + return nothing |
| 25 | +end |
| 26 | + |
| 27 | +function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::AbstractAlgorithm) |
| 28 | + m = LinearAlgebra.checksquare(A) |
| 29 | + @check_size(expA, (m, m)) |
| 30 | + @check_scalar(expA, A, (τ isa Real) ? identity : complex) |
| 31 | + return nothing |
| 32 | +end |
| 33 | + |
| 34 | +function check_input(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, ::DiagonalAlgorithm) |
| 35 | + m = LinearAlgebra.checksquare(A) |
| 36 | + @assert isdiag(A) |
| 37 | + @assert expA isa Diagonal |
| 38 | + @check_size(expA, (m, m)) |
| 39 | + @check_scalar(expA, A, (τ isa Real) ? identity : complex) |
| 40 | + return nothing |
| 41 | +end |
| 42 | + |
| 43 | +# Outputs |
| 44 | +# ------- |
| 45 | +initialize_output(::typeof(exponential!), A::AbstractMatrix, ::AbstractAlgorithm) = A |
| 46 | +initialize_output(::typeof(exponential!), (τ, A)::Tuple{T, AbstractMatrix}, ::AbstractAlgorithm) where {T <: Real} = A |
| 47 | +initialize_output(::typeof(exponential!), (τ, A)::Tuple{Number, AbstractMatrix}, ::AbstractAlgorithm) = complex(A) |
| 48 | + |
| 49 | +# Implementation |
| 50 | +# -------------- |
| 51 | +function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaLA) |
| 52 | + check_input(exponential!, A, expA, alg) |
| 53 | + A = LinearAlgebra.exp!(A) |
| 54 | + A === expA || copy!(expA, A) |
| 55 | + return expA |
| 56 | +end |
| 57 | + |
| 58 | +exponential!(A, expA, alg::MatrixFunctionViaEigh) = exponential!((one(eltype(A)), A), expA, alg) |
| 59 | +exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEig) = exponential!((one(eltype(A)), A), expA, alg) |
| 60 | + |
| 61 | +function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::AbstractAlgorithm) |
| 62 | + expA .= A .* τ |
| 63 | + return exponential!(expA, expA, alg) |
| 64 | +end |
| 65 | + |
| 66 | +function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::MatrixFunctionViaEigh) |
| 67 | + check_input(exponential!, (τ, A), expA, alg) |
| 68 | + D, V = eigh_full!(A, alg.eigh_alg) |
| 69 | + if eltype(A) <: Real |
| 70 | + if eltype(τ) <: Real |
| 71 | + VexpD = rmul!(V, exponential!((τ / 2, D), D)) |
| 72 | + else |
| 73 | + VexpD = V * exponential((τ / 2, D)) |
| 74 | + end |
| 75 | + return mul!(expA, VexpD, transpose(VexpD)) |
| 76 | + else |
| 77 | + if eltype(τ) <: Real |
| 78 | + VexpD = V * exponential!((τ, D), D) |
| 79 | + else |
| 80 | + VexpD = V * exponential((τ, D)) |
| 81 | + end |
| 82 | + return mul!(expA, VexpD, V') |
| 83 | + end |
| 84 | +end |
| 85 | + |
| 86 | +function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::MatrixFunctionViaEig) |
| 87 | + check_input(exponential!, (τ, A), expA, alg) |
| 88 | + D, V = eig_full!(A, alg.eig_alg) |
| 89 | + if eltype(A) <: Real && eltype(τ) <: Real |
| 90 | + VexpD = V * exponential!((τ, D), D) |
| 91 | + expAc = rdiv!(VexpD, LinearAlgebra.lu!(V)) |
| 92 | + return expA .= real.(expAc) |
| 93 | + else |
| 94 | + expA .= V .* transpose(diagview(exponential!((τ, D), D))) |
| 95 | + return rdiv!(expA, LinearAlgebra.lu!(V)) |
| 96 | + end |
| 97 | +end |
| 98 | + |
| 99 | +# Diagonal logic |
| 100 | +# -------------- |
| 101 | +function exponential!(A, expA, alg::DiagonalAlgorithm) |
| 102 | + check_input(exponential!, A, expA, alg) |
| 103 | + return map_diagonal!(exp, expA, A) |
| 104 | +end |
| 105 | + |
| 106 | +function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::DiagonalAlgorithm) |
| 107 | + check_input(exponential!, (τ, A), expA, alg) |
| 108 | + return map_diagonal!(x -> exp(x * τ), expA, A) |
| 109 | +end |
0 commit comments