Apply JuliaFormatter to fix code formatting

ext/ExponentialUtilitiesStaticArraysExt.jl

@@ -1,98 +1,106 @@
 module ExponentialUtilitiesStaticArraysExt
 
-export default_tolerance,theta,THETA32,THETA64
+export default_tolerance, theta, THETA32, THETA64
 
 using StaticArrays
 import Base: @propagate_inbounds
-import LinearAlgebra: tr,I,opnorm,norm
+import LinearAlgebra: tr, I, opnorm, norm
 import ExponentialUtilities
 
 # Look-Up Table Generation
 default_tolerance(::Type{T}) where {T <: AbstractFloat} = eps(T)/2
-@inline function trexp(M::Integer,x::T) where {T}
+@inline function trexp(M::Integer, x::T) where {T}
     y = T <: BigInt ? one(BigFloat) : T <: Integer ? one(Float64) : one(T)
-    for m ∈ M:-1:1
+    for m in M:-1:1
         y = 1+x/m*y
     end
     return y
 end
-h(M::Integer,x::Number) = log(exp(-x)*trexp(M,x))
-h̃(M::Integer,x::Number) = ifelse(isodd(M),-1,1)*h(M,-x)
-function θf((M,ϵ)::Tuple{<:Integer,<:Number},x::Number)
-    return h̃(M+1,x)/x-ϵ
+h(M::Integer, x::Number) = log(exp(-x)*trexp(M, x))
+h̃(M::Integer, x::Number) = ifelse(isodd(M), -1, 1)*h(M, -x)
+function θf((M, ϵ)::Tuple{<:Integer, <:Number}, x::Number)
+    return h̃(M+1, x)/x-ϵ
 end
-θf(M::Integer,ϵ::Number) = Base.Fix1(θf,(M,ϵ))
-function θfp(M::Integer,x::Number)
-    Tk = trexp(M+1,-x)
-    Tkm1 = trexp(M,-x)
-    return ifelse(isodd(M),-1,1)/x^2*(log(Tk)+x*Tkm1/Tk)
+θf(M::Integer, ϵ::Number) = Base.Fix1(θf, (M, ϵ))
+function θfp(M::Integer, x::Number)
+    Tk = trexp(M+1, -x)
+    Tkm1 = trexp(M, -x)
+    return ifelse(isodd(M), -1, 1)/x^2*(log(Tk)+x*Tkm1/Tk)
 end
-θfp(M::Integer) = Base.Fix1(θfp,M)
+θfp(M::Integer) = Base.Fix1(θfp, M)
 
-function newton_find_zero(f::Function,dfdx::Function,x0::Real;xrtol::Real=eps(typeof(x0))/2,maxiter::Integer=100)
-    0 ≤ xrtol ≤ 1 || throw(DomainError(xrtol,"relative tolerance in x must be in [0,1]"))
-    maxiter > 0 || throw(DomainError(maxiter,"maxiter should be a positive integer"))
+function newton_find_zero(f::Function, dfdx::Function, x0::Real;
+        xrtol::Real = eps(typeof(x0))/2, maxiter::Integer = 100)
+    0 ≤ xrtol ≤ 1 || throw(DomainError(xrtol, "relative tolerance in x must be in [0,1]"))
+    maxiter > 0 || throw(DomainError(maxiter, "maxiter should be a positive integer"))
     x, xp = x0, typemax(x0)
-    for _ ∈ 1:maxiter
+    for _ in 1:maxiter
         xp = x
         x -= f(x)/dfdx(x)
-        if abs(x-xp) ≤ xrtol*max(x,xp) || !isfinite(x)
+        if abs(x-xp) ≤ xrtol*max(x, xp) || !isfinite(x)
             break
         end
     end
     return x
 end
-function calc_thetas(m_max::Integer,::Type{T};tol::T=default_tolerance(T)) where {T <: AbstractFloat}
-    m_max > 0 || throw(DomainError(m_max,"argument m_max must be positive"))
+function calc_thetas(m_max::Integer, ::Type{T}; tol::T = default_tolerance(T)) where {T <:
+                                                                                      AbstractFloat}
+    m_max > 0 || throw(DomainError(m_max, "argument m_max must be positive"))
     ϵ = BigFloat(tol)
-    θ = Vector{T}(undef,m_max+1)
+    θ = Vector{T}(undef, m_max+1)
     @inbounds θ[1] = eps(T)
-    @inbounds for m=1:m_max
-        θ[m+1] = newton_find_zero(θf(m,ϵ),θfp(m),big(θ[m]),xrtol=ϵ)
+    @inbounds for m in 1:m_max
+        θ[m + 1] = newton_find_zero(θf(m, ϵ), θfp(m), big(θ[m]), xrtol = ϵ)
     end
     return θ
 end
 
 const P_MAX = 8
 const M_MAX = 55
-const THETA32 = Tuple(calc_thetas(M_MAX,Float32))
-const THETA64 = Tuple(calc_thetas(M_MAX,Float64))
+const THETA32 = Tuple(calc_thetas(M_MAX, Float32))
+const THETA64 = Tuple(calc_thetas(M_MAX, Float64))
 
-@propagate_inbounds theta(::Type{Float64},m::Integer) = THETA64[m]
-@propagate_inbounds theta(::Type{Float32},m::Integer) = THETA32[m]
-@propagate_inbounds theta(::Type{Complex{T}},m::Integer) where {T} = theta(T,m)
-@propagate_inbounds theta(::Type{T},::Integer) where {T} = throw(DomainError(T,"type must be either Float32 or Float64"))
-@propagate_inbounds theta(x::Number,m::Integer) = theta(typeof(x),m)
+@propagate_inbounds theta(::Type{Float64}, m::Integer) = THETA64[m]
+@propagate_inbounds theta(::Type{Float32}, m::Integer) = THETA32[m]
+@propagate_inbounds theta(::Type{Complex{T}}, m::Integer) where {T} = theta(T, m)
+@propagate_inbounds theta(
+    ::Type{T}, ::Integer) where {T} = throw(DomainError(T, "type must be either Float32 or Float64"))
+@propagate_inbounds theta(x::Number, m::Integer) = theta(typeof(x), m)
 
 # runtime parameter search
-@propagate_inbounds @inline function calculate_s(α::T,m::I)::I where {T <: Number,I <: Integer}
-    return ceil(I,α/theta(T,m))
+@propagate_inbounds @inline function calculate_s(α::T, m::I)::I where {
+        T <: Number, I <: Integer}
+    return ceil(I, α/theta(T, m))
 end
-@propagate_inbounds @inline function parameter_search(nA::Number,m::I)::I where {I <: Integer}
-    return m*calculate_s(nA,m)
+@propagate_inbounds @inline function parameter_search(nA::Number, m::I)::I where {I <:
+                                                                                  Integer}
+    return m*calculate_s(nA, m)
 end
-@propagate_inbounds @inline function parameters(A::SMatrix{N,N,T})::Tuple{Int,Int} where {N,T}
-    1 ≤ N ≤ 50 || throw(DomainError(N,"leading dimension of A must be ≤ 50; larger matrices require Higham's 1-norm estimation algorithm"))
-    nA = opnorm(A,1)
-    iszero(nA) && return (0,1)
-    @inbounds if nA ≤ 4theta(T,M_MAX)*P_MAX*(P_MAX+3)/(M_MAX*1)
-        mo = argmin(Base.Fix1(parameter_search,nA),1:M_MAX)
-        s = calculate_s(nA,mo)
-        return (mo,s)
+@propagate_inbounds @inline function parameters(A::SMatrix{
+        N, N, T})::Tuple{Int, Int} where {N, T}
+    1 ≤ N ≤ 50 || throw(DomainError(N,
+        "leading dimension of A must be ≤ 50; larger matrices require Higham's 1-norm estimation algorithm"))
+    nA = opnorm(A, 1)
+    iszero(nA) && return (0, 1)
+    @inbounds if nA ≤ 4theta(T, M_MAX)*P_MAX*(P_MAX+3)/(M_MAX*1)
+        mo = argmin(Base.Fix1(parameter_search, nA), 1:M_MAX)
+        s = calculate_s(nA, mo)
+        return (mo, s)
     else
         Aᵐ = A*A
-        pη = √(opnorm(Aᵐ,1))
-        (Cmo::Int,mo::Int) = (typemax(Int),1)
-        for p ∈ 2:P_MAX
+        pη = √(opnorm(Aᵐ, 1))
+        (Cmo::Int, mo::Int) = (typemax(Int), 1)
+        for p in 2:P_MAX
             Aᵐ *= A
-            η = opnorm(Aᵐ,1)^inv(p+1)
-            α = max(pη,η)
+            η = opnorm(Aᵐ, 1)^inv(p+1)
+            α = max(pη, η)
             pη = η
-            (Cmp::Int,mp::Int) = findmin(Base.Fix1(parameter_search,α),p*(p-1)-1:M_MAX)
-            (Cmo,mo) = min((Cmp,mp),(Cmo,mo))
+            (Cmp::Int,
+                mp::Int) = findmin(Base.Fix1(parameter_search, α), (p * (p - 1) - 1):M_MAX)
+            (Cmo, mo) = min((Cmp, mp), (Cmo, mo))
         end
-        s = max(Cmo÷mo,1)
-        return (mo,s)
+        s = max(Cmo÷mo, 1)
+        return (mo, s)
     end
 end
 
@@ -105,11 +113,12 @@ end
     Presently, the relative tolerance is fixed to eps(T)/2 where T is the type of the
     output.
 """
-@propagate_inbounds function ExponentialUtilities.expv(t::Number,A::SMatrix{N,N,T},v::SVector{N}; kwarg...) where {N,T}
-    Ti = promote_type(StaticArrays.arithmetic_closure(T),eltype(v))
+@propagate_inbounds function ExponentialUtilities.expv(
+        t::Number, A::SMatrix{N, N, T}, v::SVector{N}; kwarg...) where {N, T}
+    Ti = promote_type(StaticArrays.arithmetic_closure(T), eltype(v))
     N ≤ 4 && return exp(t*A)*v
-    Ai::SMatrix{N,N,Ti} = A
-    
+    Ai::SMatrix{N, N, Ti} = A
+
     μ = tr(Ai)/N
     Ai -= μ*I
     Ai *= t
@@ -118,18 +127,18 @@ end
     Ai /= s
     η = exp(μ*t/s)
     ϵ = default_tolerance(T)
-    for _ ∈ 1:s
-        c₁ = norm(v,Inf)
-        for j ∈ 1:mo
+    for _ in 1:s
+        c₁ = norm(v, Inf)
+        for j in 1:mo
             v = (Ai*v)/j
             F += v
-            c₂ = norm(v,Inf)
-            c₁+c₂ ≤ ϵ*norm(F,Inf) && break
+            c₂ = norm(v, Inf)
+            c₁+c₂ ≤ ϵ*norm(F, Inf) && break
             c₁ = c₂
         end
         F *= η
         v = F
-        all(isfinite,v) || break
+        all(isfinite, v) || break
     end
 
     return F

src/arnoldi.jl

@@ -12,11 +12,10 @@ The dimension of the subspace, `Ks.m`, can be dynamically altered but should
 be smaller than `maxiter`, the maximum allowed arnoldi iterations.
 
 The type of the (extended) orthonormal basis vector matrix `V` may be specified
-as `VType`. This is required e. g. for `GPUArray`s. 
+as `VType`. This is required e. g. for `GPUArray`s.
 
 `U` determines `eltype(H)`.
 
-
     getV(Ks) -> V
     getH(Ks) -> H
 
@@ -50,7 +49,7 @@ function KrylovSubspace{T, U, VType}(n::Integer, maxiter::Integer = 30,
 end
 
 KrylovSubspace{T}(args...) where {T} = KrylovSubspace{T, T}(args...)
-KrylovSubspace{T,U}(args...) where {T,U} = KrylovSubspace{T, U, Matrix{T}}(args...)
+KrylovSubspace{T, U}(args...) where {T, U} = KrylovSubspace{T, U, Matrix{T}}(args...)
 
 getV(Ks::KrylovSubspace) = @view(Ks.V[:, 1:(Ks.m + 1)])
 getH(Ks::KrylovSubspace) = @view(Ks.H[1:(Ks.m + 1), 1:(Ks.m + !iszero(Ks.augmented))])

src/exp_baseexp.jl

@@ -45,14 +45,14 @@ function exponential!(A::StridedMatrix{T}, method::ExpMethodHigham2005Base,
     if (nA <= 2.1)
         if nA > 0.95
             C = T[17643225600.0, 8821612800.0, 2075673600.0, 302702400.0,
-                30270240.0, 2162160.0, 110880.0, 3960.0,
-                90.0, 1.0]
+            30270240.0, 2162160.0, 110880.0, 3960.0,
+            90.0, 1.0]
         elseif nA > 0.25
             C = T[17297280.0, 8648640.0, 1995840.0, 277200.0,
-                25200.0, 1512.0, 56.0, 1.0]
+            25200.0, 1512.0, 56.0, 1.0]
         elseif nA > 0.015
             C = T[30240.0, 15120.0, 3360.0,
-                420.0, 30.0, 1.0]
+            420.0, 30.0, 1.0]
         else
             C = T[120.0, 60.0, 12.0, 1.0]
         end
@@ -78,10 +78,10 @@ function exponential!(A::StridedMatrix{T}, method::ExpMethodHigham2005Base,
             A ./= convert(T, 2^si)
         end
         C = T[64764752532480000.0, 32382376266240000.0, 7771770303897600.0,
-            1187353796428800.0, 129060195264000.0, 10559470521600.0,
-            670442572800.0, 33522128640.0, 1323241920.0,
-            40840800.0, 960960.0, 16380.0,
-            182.0, 1.0]
+        1187353796428800.0, 129060195264000.0, 10559470521600.0,
+        670442572800.0, 33522128640.0, 1323241920.0,
+        40840800.0, 960960.0, 16380.0,
+        182.0, 1.0]
         mul!(A2, A, A)
         @. U = C[2] * P
         @. V = C[1] * P

src/exp_generic.jl

@@ -97,7 +97,7 @@ _const(A::Array) = Base.Experimental.Const(A)
                         for k in Kaxis
                             Base.Cartesian.@nexprs $NU j->Cmn_j = muladd(_const(A)[mm, k],
                                 _const(B)[k,
-                                    n + j],
+                                n + j],
                                 Cmn_j)
                         end
                         Base.Cartesian.@nexprs $NU j->C[mm, n + j] = Cmn_j

src/exp_sparse.jl

@@ -1,9 +1,10 @@
 
-for expmeth in [ExpMethodDiagonalization, ExpMethodGeneric, ExpMethodHigham2005, ExpMethodHigham2005Base, ExpMethodNative]    
-    @eval function exponential!(A::AbstractSparseArray, method::$expmeth, cache=nothing)
-        throw(ErrorException("exp(A) on a sparse matrix is generally dense. This operation is "*
-        "not allowed with exponential. If you wished to compute exp(At)*v, see expv. "*
-        "Otherwise to override this error, densify the matrix before calling, "*
-        "i.e. exponential!(Matrix(A))"))
+for expmeth in [ExpMethodDiagonalization, ExpMethodGeneric,
+    ExpMethodHigham2005, ExpMethodHigham2005Base, ExpMethodNative]
+    @eval function exponential!(A::AbstractSparseArray, method::$expmeth, cache = nothing)
+        throw(ErrorException("exp(A) on a sparse matrix is generally dense. This operation is " *
+                             "not allowed with exponential. If you wished to compute exp(At)*v, see expv. " *
+                             "Otherwise to override this error, densify the matrix before calling, " *
+                             "i.e. exponential!(Matrix(A))"))
     end
-end
+end

src/kiops.jl

@@ -210,7 +210,11 @@ function kiops(tau_out, A, u; mmin::Int = 10, mmax::Int = 128, m::Int = min(mmin
 
         # Check error against target
         if omega <= delta
-            tau_now, l, j, reject, ireject, step = kiops_update_solution!(tau_now, tau,
+            tau_now, l,
+            j,
+            reject,
+            ireject,
+            step = kiops_update_solution!(tau_now, tau,
                 tau_out, w, l, V,
                 F, H, beta, j, n,
                 step, numSteps,

src/krylov_phiv.jl

@@ -12,7 +12,7 @@ function Base.resize!(C::ExpvCache{T}, maxiter::Int) where {T}
 end
 function get_cache(C::ExpvCache, m::Int)
     m^2 > length(C.mem) && resize!(C, m) # resize the cache if needed
-    reshape(@view(C.mem[1:(m^2)]), m, m)
+    reshape(@view(C.mem[1:(m ^ 2)]), m, m)
 end
 
 ############################
@@ -195,15 +195,15 @@ function get_caches(C::PhivCache{useview, T}, m::Int, p::Int) where {useview, T}
     offset = m
 
     if useview
-        Hcopy = reshape(@view(C.mem[(offset + 1):(offset + m^2)]), m, m)
+        Hcopy = reshape(@view(C.mem[(offset + 1):(offset + m ^ 2)]), m, m)
         offset += m^2
-        C1 = reshape(@view(C.mem[(offset + 1):(offset + (m + p)^2)]), m + p, m + p)
+        C1 = reshape(@view(C.mem[(offset + 1):(offset + (m + p) ^ 2)]), m + p, m + p)
         offset += (m + p)^2
         C2 = reshape(@view(C.mem[(offset + 1):(offset + m * (p + 1))]), m, p + 1)
     else
-        Hcopy = reshape(C.mem[(offset + 1):(offset + m^2)], m, m)
+        Hcopy = reshape(C.mem[(offset + 1):(offset + m ^ 2)], m, m)
         offset += m^2
-        C1 = reshape(C.mem[(offset + 1):(offset + (m + p)^2)], m + p, m + p)
+        C1 = reshape(C.mem[(offset + 1):(offset + (m + p) ^ 2)], m + p, m + p)
         offset += (m + p)^2
         C2 = reshape(C.mem[(offset + 1):(offset + m * (p + 1))], m, p + 1)
     end

src/krylov_phiv_adaptive.jl

@@ -183,7 +183,8 @@ function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractM
         if Ks.wasbreakdown
             tau = tend - t
         end
-        _, epsilon = phiv!(P, tau, Ks, p + 1; cache = phiv_cache, correct = correct,
+        _,
+        epsilon = phiv!(P, tau, Ks, p + 1; cache = phiv_cache, correct = correct,
             errest = true)
         verbose && println("t = $t, m = $m, tau = $tau, error estimate = $epsilon")
         if adaptive
@@ -195,7 +196,9 @@ function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractM
             kappa = 2.0
             maxtau = tend - t
             while omega > delta # inner loop of Algorithm 3
-                m_new, tau_new, q, kappa = _phiv_timestep_adapt(m, tau, epsilon, m_old,
+                m_new, tau_new,
+                q,
+                kappa = _phiv_timestep_adapt(m, tau, epsilon, m_old,
                     tau_old, epsilon_old, q,
                     kappa,
                     gamma, omega, maxtau, n, p,
@@ -208,7 +211,8 @@ function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractM
                 # Compute ϕp(tau*A)wp using the new parameters
                 arnoldi!(Ks, A, @view(W[:, end]); tol = tol, m = m, opnorm = opnorm,
                     iop = iop)
-                _, epsilon_new = phiv!(P, tau, Ks, p + 1; cache = phiv_cache,
+                _,
+                epsilon_new = phiv!(P, tau, Ks, p + 1; cache = phiv_cache,
                     correct = correct, errest = true)
                 epsilon, epsilon_old = epsilon_new, epsilon
                 omega = (tend / tau) * (epsilon / abstol)

src/phi.jl

@@ -143,7 +143,7 @@ function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches =
 end
 function phi!(out::Vector{Diagonal{T, V}}, A::Diagonal{T, V}, k::Integer;
         caches = nothing) where {T <: Number, V <: AbstractVector{T}}
-    for i in axes(A,1)
+    for i in axes(A, 1)
         phiz = phi(A[i, i], k; cache = caches)
         for j in 1:(k + 1)
             out[j][i, i] = phiz[j]

test/basictests.jl

@@ -126,7 +126,7 @@ end
 end
 
 @testset "naive_matmul" begin
-    A = Matrix(reshape((1.0:(23.0^2)) ./ 700, (23, 23)))
+    A = Matrix(reshape((1.0:(23.0 ^ 2)) ./ 700, (23, 23)))
     @test exp_generic(A) ≈ exp(A)
     # FiniteDiff.finite_difference_gradient(sum ∘ exp, A)
     @test ForwardDiff.gradient(sum ∘ exp_generic, A) ≈
@@ -215,10 +215,11 @@ end
 
 @testset "Static Arrays" begin
     Random.seed!(0)
-    for N in (3,4,6,8),t in (0.1,1.0,10.0)
-        A = I+randn(SMatrix{N,N,Float64})/3
-        b = randn(SVector{N,Float64})
-        @test expv(t,A,b) ≈ exp(t*A)*b
+    for N in (3, 4, 6, 8), t in (0.1, 1.0, 10.0)
+
+        A = I+randn(SMatrix{N, N, Float64})/3
+        b = randn(SVector{N, Float64})
+        @test expv(t, A, b) ≈ exp(t*A)*b
     end
 end
 
@@ -293,6 +294,7 @@ end
         rand(n, n)
     ]
         for b in [rand(ComplexF64, n), rand(n)], t in [1e-2, 1e-2im, 1e-2 + 1e-2im]
+
             @test exp(t * A) * b ≈ expv(t, A, b; m = m)
         end
     end