Run formatter on entire repo (#123)

bench/example.jl

@@ -9,8 +9,8 @@ using Arpack
 
 function helloworld(n = 40, nev = 4)
     A = randn(n, n)
-    @time S, hist = partialschur(A, nev=nev, which = LR(), tol = 1e-6)
-    @time eigs(A, nev=nev, which = :LR, tol = 1e-6)
+    @time S, hist = partialschur(A, nev = nev, which = LR(), tol = 1e-6)
+    @time eigs(A, nev = nev, which = :LR, tol = 1e-6)
     @show norm(A * S.Q - S.Q * S.R)
     println(hist)
     θs = eigvals(S.R)

bench/partial_schur.jl

@@ -6,39 +6,47 @@ using LinearAlgebra
 using SparseArrays
 using DelimitedFiles
 
-mymatrix(n) = spdiagm(-1 => fill(-1.0, n-1), 0 => fill(2.0, n), 1 => fill(-1.001, n-1))
+mymatrix(n) = spdiagm(-1 => fill(-1.0, n - 1), 0 => fill(2.0, n), 1 => fill(-1.001, n - 1))
 
 struct ShiftAndInvert{TA,TB,TT}
     A_lu::TA
     B::TB
     temp::TT
 end
 
-function (M::ShiftAndInvert)(y,x)
+function (M::ShiftAndInvert)(y, x)
     mul!(M.temp, M.B, x)
     # @show typeof(y), typeof(M.A_lu), typeof(M.temp)
     ldiv!(y, M.A_lu, M.temp)
 end
 
-function construct_linear_map(A,B)
-    a = ShiftAndInvert(factorize(A),B,Vector{eltype(A)}(undef, size(A,1)))
-    LinearMap(a, size(A,1), ismutating=true)
+function construct_linear_map(A, B)
+    a = ShiftAndInvert(factorize(A), B, Vector{eltype(A)}(undef, size(A, 1)))
+    LinearMap(a, size(A, 1), ismutating = true)
 end
 
 # For inverting in non-generealized case
 function construct_linear_map(A)
     a = factorize(A)
-    function f(y,x)
+    function f(y, x)
         ldiv!(y, a, x)
     end
-    LinearMap{eltype(A)}(f, size(A,1))
+    LinearMap{eltype(A)}(f, size(A, 1))
 end
 
 function bencharnoldi()
     A = mymatrix(6000)
     A = construct_linear_map(mymatrix(6000))
     target = LM()
 
-    @benchmark partialschur(B, mindim=11, maxdim=22, nev=10, tol=1e-10, restarts=100000, which=tar) setup = (B = $A; tar = $target)
+    @benchmark partialschur(
+        B,
+        mindim = 11,
+        maxdim = 22,
+        nev = 10,
+        tol = 1e-10,
+        restarts = 100000,
+        which = tar,
+    ) setup = (B = $A; tar = $target)
 
 end

bench/schur.jl

@@ -4,7 +4,8 @@ using BenchmarkTools
 function bench_schur(n = 30)
     H = triu(rand(n, n), -1)
     Q = Matrix{Float64}(I, n, n)
-    a = @benchmark ArnoldiMethod.local_schurfact!(HH, 1, $n, QQ) setup = (HH = copy($H); QQ = copy($Q))
+    a = @benchmark ArnoldiMethod.local_schurfact!(HH, 1, $n, QQ) setup =
+        (HH = copy($H); QQ = copy($Q))
     b = @benchmark eigvals(HH) setup = (HH = copy($H))
     a, b
 end

docs/make.jl

@@ -5,12 +5,12 @@ makedocs(
     sitename = "ArnoldiMethod.jl",
     format = Documenter.HTML(),
     pages = [
-	"Home" => "index.md",
-	"Theory" => "theory.md",
-	"Using ArnoldiMethod.jl" => [
-	    "Getting started" => "usage/01_getting_started.md",
-	    "Transformations" => "usage/02_spectral_transformations.md"
-	]
+        "Home" => "index.md",
+        "Theory" => "theory.md",
+        "Using ArnoldiMethod.jl" => [
+            "Getting started" => "usage/01_getting_started.md",
+            "Transformations" => "usage/02_spectral_transformations.md",
+        ],
     ],
     warnonly = [:missing_docs, :cross_references],
 )
@@ -20,5 +20,5 @@ makedocs(
 # for more information.
 deploydocs(
     repo = "github.com/JuliaLinearAlgebra/ArnoldiMethod.jl.git",
-    devbranch="master"
+    devbranch = "master",
 )

src/ArnoldiMethod.jl

@@ -21,9 +21,10 @@ struct Arnoldi{T,TV<:StridedMatrix{T},TH<:StridedMatrix{T}}
     H::TH
 
     function Arnoldi{T}(matrix_order::Int, krylov_dimension::Int) where {T}
-        krylov_dimension <= matrix_order || throw(ArgumentError("Krylov dimension should be less than matrix order."))
+        krylov_dimension <= matrix_order ||
+            throw(ArgumentError("Krylov dimension should be less than matrix order."))
         V = Matrix{T}(undef, matrix_order, krylov_dimension + 1)
-        H = zeros(T, krylov_dimension + 1, krylov_dimension)    
+        H = zeros(T, krylov_dimension + 1, krylov_dimension)
         return new{T,typeof(V),typeof(H)}(V, H)
     end
 end
@@ -83,4 +84,4 @@ include("eigenvector_uppertriangular.jl")
 include("show.jl")
 
 
-end
+end

src/eigenvector_uppertriangular.jl

@@ -6,51 +6,51 @@ Solve the problem (R[1:k,1:k] - λI) \\ x[1:k] in-place.
 function shifted_backward_sub!(x, R::AbstractMatrix{Tm}, λ, k) where {Tm<:Real}
     # Real arithmetic with quasi-upper triangular R
     @inbounds while k > 0
-        if k > 1 && R[k,k-1] != zero(Tm)
+        if k > 1 && R[k, k-1] != zero(Tm)
             # Solve 2x2 problem
-            R11, R12 = R[k-1,k-1] - λ, R[k-1,k]
-            R21, R22 = R[k  ,k-1]    , R[k  ,k] - λ
+            R11, R12 = R[k-1, k-1] - λ, R[k-1, k]
+            R21, R22 = R[k, k-1], R[k, k] - λ
             det = R11 * R22 - R21 * R12
-            a1 = ( R22 * x[k-1] - R12 * x[k-0]) / det
+            a1 = (R22 * x[k-1] - R12 * x[k-0]) / det
             a2 = (-R21 * x[k-1] + R11 * x[k-0]) / det
             x[k-1] = a1
             x[k-0] = a2
 
             # Backward substitute
-            for i = 1 : k-2
-                x[i] -= R[i,k-1] * x[k-1] + R[i,k-0] * x[k-0]
+            for i = 1:k-2
+                x[i] -= R[i, k-1] * x[k-1] + R[i, k-0] * x[k-0]
             end
             k -= 2
         else
             # Solve 1x1 "problem"
-            x[k] /= R[k,k] - λ
-    
+            x[k] /= R[k, k] - λ
+
             # Backward substitute
-            for i = 1 : k - 1
-                x[i] -= R[i,k] * x[k]
+            for i = 1:k-1
+                x[i] -= R[i, k] * x[k]
             end
 
             k -= 1
         end
     end
-    
+
     x
 end
 
 function shifted_backward_sub!(x, R::AbstractMatrix, λ, k)
     # Generic implementation, upper triangular R
     @inbounds while k > 0
         # Solve 1x1 "problem"
-        x[k] /= R[k,k] - λ
+        x[k] /= R[k, k] - λ
 
         # Backward substitute
-        for i = 1 : k - 1
-            x[i] -= R[i,k] * x[k]
+        for i = 1:k-1
+            x[i] -= R[i, k] * x[k]
         end
 
         k -= 1
     end
-    
+
     x
 end
 
@@ -60,40 +60,44 @@ end
 Store the `j`th eigenvector of an upper triangular matrix `R` in `x`.
 In the end `norm(x[1:k]) = 1` This function leaves x[k+1:end] untouched!
 """
-function collect_eigen!(x::AbstractVector{Tv}, R::AbstractMatrix{Tm}, j::Integer) where {Tm<:Real,Tv<:Number}
+function collect_eigen!(
+    x::AbstractVector{Tv},
+    R::AbstractMatrix{Tm},
+    j::Integer,
+) where {Tm<:Real,Tv<:Number}
     n = size(R, 2)
 
     @inbounds begin
         # If it's a conjugate pair with the next index, then just increment j.
-        if j < n && R[j+1,j] != zero(Tv)
+        if j < n && R[j+1, j] != zero(Tv)
             j += 1
         end
 
         # Initialize the rhs and do the first backward substitution
         # Then do the rest of the shifted backward substitution in another function,
         # cause λ is either real or complex.
-        if j > 1 && R[j,j-1] != 0
+        if j > 1 && R[j, j-1] != 0
             # Complex arithmetic
-            R11, R21 = R[j-1,j-1], R[j-0,j-1]
-            R12, R22 = R[j-1,j-0], R[j-0,j-0]
+            R11, R21 = R[j-1, j-1], R[j-0, j-1]
+            R12, R22 = R[j-1, j-0], R[j-0, j-0]
             det = R11 * R22 - R21 * R12
             tr = R11 + R22
             λ = (tr + sqrt(complex(tr * tr - 4det))) / 2
             x[j-1] = -R12 / (R11 - λ)
             x[j-0] = one(Tv)
-            for i = 1 : j-2
-                x[i] = -R[i,j-1] * x[j-1] - R[i,j]
+            for i = 1:j-2
+                x[i] = -R[i, j-1] * x[j-1] - R[i, j]
             end
 
-            shifted_backward_sub!(x, R, λ, j-2)
+            shifted_backward_sub!(x, R, λ, j - 2)
         else
             # Real arithmetic
-            λ = R[j,j]
+            λ = R[j, j]
             x[j] = one(Tv)
-            for i = 1 : j-1
-                x[i] = -R[i,j]
+            for i = 1:j-1
+                x[i] = -R[i, j]
             end
-            shifted_backward_sub!(x, R, λ, j-1)
+            shifted_backward_sub!(x, R, λ, j - 1)
         end
 
         # Normalize
@@ -114,13 +118,13 @@ function collect_eigen!(x::AbstractVector{Tv}, R::AbstractMatrix, j::Integer) wh
     n = size(R, 2)
 
     @inbounds begin
-        λ = R[j,j]
+        λ = R[j, j]
         x[j] = one(Tv)
-        for i = 1 : j-1
-            x[i] = -R[i,j]
+        for i = 1:j-1
+            x[i] = -R[i, j]
         end
 
-        shifted_backward_sub!(x, R, λ, j-1)
+        shifted_backward_sub!(x, R, λ, j - 1)
 
         # Normalize
         nrm = zero(real(Tv))
@@ -134,4 +138,4 @@ function collect_eigen!(x::AbstractVector{Tv}, R::AbstractMatrix, j::Integer) wh
     end
 
     return j
-end
+end

src/eigvals.jl

@@ -3,7 +3,12 @@
 
 Puts the eigenvalues of a quasi-upper triangular matrix A in the λs vector.
 """
-function copy_eigenvalues!(λs, A::AbstractMatrix{T}, range = OneTo(size(A, 2)), tol = eps(real(T))) where {T}
+function copy_eigenvalues!(
+    λs,
+    A::AbstractMatrix{T},
+    range = OneTo(size(A, 2)),
+    tol = eps(real(T)),
+) where {T}
     i = first(range)
 
     @inbounds while i < last(range)
@@ -12,17 +17,17 @@ function copy_eigenvalues!(λs, A::AbstractMatrix{T}, range = OneTo(size(A, 2)),
             i += 1
         else
             # Conjugate pair
-            d = A[i,i] * A[i+1,i+1] - A[i,i+1] * A[i+1,i]
-            x = (A[i,i] + A[i+1,i+1]) / 2
-            y = sqrt(complex(x*x - d))
-            λs[i + 0] = x + y
-            λs[i + 1] = x - y
+            d = A[i, i] * A[i+1, i+1] - A[i, i+1] * A[i+1, i]
+            x = (A[i, i] + A[i+1, i+1]) / 2
+            y = sqrt(complex(x * x - d))
+            λs[i+0] = x + y
+            λs[i+1] = x - y
             i += 2
         end
     end
 
     @inbounds if i == last(range)
-        λs[i] = A[i, i] 
+        λs[i] = A[i, i]
     end
 
     return λs
@@ -38,12 +43,12 @@ function eigenvalue(R, i)
     n = minimum(size(R))
 
     @inbounds begin
-        if i == n || iszero(R[i+1,i])
-            return complex(R[i,i])
+        if i == n || iszero(R[i+1, i])
+            return complex(R[i, i])
         else
-            d = R[i,i] * R[i+1,i+1] - R[i,i+1] * R[i+1,i]
-            x = (R[i,i] + R[i+1,i+1]) / 2
-            y = sqrt(complex(x*x - d))
+            d = R[i, i] * R[i+1, i+1] - R[i, i+1] * R[i+1, i]
+            x = (R[i, i] + R[i+1, i+1]) / 2
+            y = sqrt(complex(x * x - d))
             return x + y
         end
     end
@@ -82,5 +87,5 @@ upper triangular matrix `R` from the partial Schur decomposition.
 """
 function partialeigen(P::PartialSchur)
     vals, vecs = eigen(P.R)
-    return vals, P.Q*vecs
+    return vals, P.Q * vecs
 end