Merge pull request #105 from fgerick/master

.github/workflows/ci.yaml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - "1.0"
+          - "1.3"
           - "1"
           - "nightly"
         os:

.gitignore

@@ -1,2 +1,3 @@
 .DS_Store
-docs/build
+docs/build
+Manifest.toml

Project.toml

@@ -1,7 +1,7 @@
 name = "ArnoldiMethod"
 uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
 author = ["Harmen Stoppels <[email protected]>", "Lauri Nyman"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -10,11 +10,12 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 StaticArrays = "0.8, 0.9, 0.10, 0.11, 0.12, 1.0"
-julia = "1"
+julia = "1.3"
 
 [extras]
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SparseArrays"]
+test = ["Test", "SparseArrays", "GenericSchur"]

src/expansion.jl

@@ -1,6 +1,5 @@
 using Random
 using LinearAlgebra
-using LinearAlgebra.BLAS: gemv!
 
 """
     reinitialize!(a::Arnoldi, j::Int = 0) → a
@@ -32,7 +31,7 @@ function reinitialize!(arnoldi::Arnoldi{T}, j::Int = 0) where {T}
 
     # Orthogonalize: h = Vprev' * v, v ← v - Vprev * Vprev' * v = v - Vprev * h
     h = Vprev' * v
-    gemv!('N', -one(T), Vprev, h, one(T), v)
+    mul!(v,Vprev,h,-one(T),one(T))
 
     # Norm after orthogonalization
     wnorm = norm(v)
@@ -41,7 +40,7 @@ function reinitialize!(arnoldi::Arnoldi{T}, j::Int = 0) where {T}
     if wnorm < η * rnorm
         rnorm = wnorm
         mul!(h, Vprev', v)
-        gemv!('N', -one(T), Vprev, h, one(T), v)
+        mul!(v,Vprev,h,-one(T),one(T))
         wnorm = norm(v)
     end
 
@@ -79,7 +78,7 @@ function orthogonalize!(arnoldi::Arnoldi{T}, j::Integer) where {T}
 
     # Orthogonalize: h = Vprev' * v, v ← v - Vprev * Vprev' * v = v - Vprev * h
     mul!(h, Vprev', v)
-    gemv!('N', -one(T), Vprev, h, one(T), v)
+    mul!(v,Vprev,h,-one(T),one(T))
 
     # Norm after orthogonalization
     wnorm = norm(v)
@@ -88,7 +87,7 @@ function orthogonalize!(arnoldi::Arnoldi{T}, j::Integer) where {T}
     if wnorm < η * rnorm
         rnorm = wnorm
         correction = Vprev' * v
-        gemv!('N', -one(T), Vprev, correction, one(T), v)
+        mul!(v,Vprev,correction,-one(T),one(T))
         h .+= correction
         wnorm = norm(v)
     end
@@ -112,7 +111,7 @@ Perform Arnoldi from `from` to `to`.
 """
 function iterate_arnoldi!(A, arnoldi::Arnoldi{T}, range::UnitRange{Int}) where {T}
     V, H = arnoldi.V, arnoldi.H
-    
+
     for j = range
         # Generate a new column of the Krylov subspace
         mul!(view(V, :, j+1), A, view(V, :,j))
@@ -127,4 +126,4 @@ function iterate_arnoldi!(A, arnoldi::Arnoldi{T}, range::UnitRange{Int}) where {
     end
 
     return arnoldi
-end
+end

src/schursort.jl

@@ -74,9 +74,13 @@ Computes the LU factorization of A using complete pivoting.
 """
 function lu(A::SMatrix{N,N,T}, ::Type{CompletePivoting}) where {N,T}
     # A, p and q should allocate, but escape analysis will eliminate this!
-    A = MMatrix(A)
-    p = @MVector fill(N, N)
+    if isbitstype(T)
+        A = MMatrix(A)
+    else
+        A = SizedMatrix(A)
+    end
     q = @MVector fill(N, N)
+    p = @MVector fill(N, N)
     singular = false
 
     # Maybe I should consider doing this recursively.
@@ -132,7 +136,11 @@ function lu(A::SMatrix{N,N,T}, ::Type{CompletePivoting}) where {N,T}
 end
 
 function (\)(LU::CompletelyPivotedLU{T,N}, b::SVector{N}) where {T,N}
-    x = MVector(b)
+    if isbitstype(T)
+        x = MVector(b)
+    else
+        x = SizedVector(b)
+    end
 
     # x ← L \ (P * b)
     for i = OneTo(N)

test/collect_eigen.jl

@@ -7,7 +7,7 @@ using LinearAlgebra
     n = 20
 
     # Exactly upper triangular matrix R
-    for T in (Float64, ComplexF64)
+    for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
 
         # Upper triangular matrix
         R = triu(rand(T, n, n))
@@ -38,37 +38,41 @@ rot(θ) = [cos(θ) sin(θ); -sin(θ) cos(θ)]
 
 @testset "Eigenvectors quasi-upper triangular matrix" begin
     # Real arithmetic with conjugate eigvals -- quasi-upper triangular R
-    n = 20
-    R = triu(rand(n, n))
-    R[1:2,1:2] .= rot(1.0) + I
-    R[10:11,10:11] .= rot(1.2) + 2I
-
-    # Compute exact eigenvectors according to LAPACK
-    if VERSION >= v"1.2.0-DEV.275"
-        λs, xs = eigen(R, sortby=nothing)
-    else
-        λs, xs = eigen(R)
-    end
+    for T in (Float64, BigFloat)
+        n = 20
+        R = triu(rand(T, n, n))
+        R[1:2,1:2] .= rot(one(T)) + I
+        R[10:11,10:11] .= rot(T(6//5)) + 2I
+
+        # Compute exact eigenvectors according to LAPACK
+        if VERSION >= v"1.2.0-DEV.275"
+            λs, xs = eigen(R, sortby=nothing)
+        else
+            λs, xs = eigen(R)
+        end
 
-    # Allocate a complex-valued eigenvector
-    x = zeros(ComplexF64, n)
+        # Allocate a complex-valued eigenvector
+        x = zeros(Complex{T}, n)
 
-    for i = 1:n
-        fill!(x, 0.0)
-        collect_eigen!(x, R, i)
+        for i = 1:n
+            fill!(x, 0)
+            collect_eigen!(x, R, i)
 
-        @test norm(x) ≈ 1
-        @test abs.(x) ≈ abs.(xs[:, i])
+            @test norm(x) ≈ 1
+            @test abs.(x) ≈ abs.(xs[:, i])
+        end
     end
 end
 
 @testset "Extract partial" begin
     n = 20
-    R = triu(rand(n, n))
-    R[1:2,1:2] .= rot(1.0) + I
-    for range = (1:3, 1:4)
-        λs = eigvals(R[range, range])
-        θs = copy_eigenvalues!(rand(ComplexF64, length(range)), R, range)
-        @test sort!(λs, by = reim) ≈ sort!(θs, by = reim)
+    for T in (Float64, BigFloat)
+        R = triu(rand(T,n, n))
+        R[1:2,1:2] .= rot(1.0) + I
+        for range = (1:3, 1:4)
+            λs = eigvals(R[range, range])
+            θs = copy_eigenvalues!(rand(Complex{T}, length(range)), R, range)
+            @test sort!(λs, by = reim) ≈ sort!(θs, by = reim)
+        end
     end
-end
+end

test/expansion.jl

@@ -12,38 +12,42 @@ end
 @testset "Arnoldi Factorization" begin
     n = 10
     max = 6
-    A = sprand(n, n, .1) + I
-
-    arnoldi = Arnoldi{Float64}(n, max)
-    reinitialize!(arnoldi)
-    V, H = arnoldi.V, arnoldi.H
-
-    # Do a few iterations
-    iterate_arnoldi!(A, arnoldi, 1:3)
-    @test A * V[:,1:3] ≈ V[:,1:4] * H[1:4,1:3]
-    @test norm(V[:,1:4]' * V[:,1:4] - I) < 1e-10
-
-    # Do the rest of the iterations.
-    iterate_arnoldi!(A, arnoldi, 4:max)
-    @test A * V[:,1:max] ≈ V * H
-    @test norm(V' * V - I) < 1e-10
+    for T in (Float64, BigFloat)
+        A = sprand(T, n, n, .1) + I
+
+        arnoldi = Arnoldi{T}(n, max)
+        reinitialize!(arnoldi)
+        V, H = arnoldi.V, arnoldi.H
+
+        # Do a few iterations
+        iterate_arnoldi!(A, arnoldi, 1:3)
+        @test A * V[:,1:3] ≈ V[:,1:4] * H[1:4,1:3]
+        @test norm(V[:,1:4]' * V[:,1:4] - I) < sqrt(eps(T))/100
+
+        # Do the rest of the iterations.
+        iterate_arnoldi!(A, arnoldi, 4:max)
+        @test A * V[:,1:max] ≈ V * H
+        @test norm(V' * V - I) < sqrt(eps(T))/100
+    end
 end
 
 @testset "Invariant subspace" begin
     # Generate a block-diagonal matrix A
-    A = [rand(4,4)  zeros(4,4); 
-         zeros(4,4) rand(4,4)]
-
-    # and an initial vector [1; 0; ... 0]
-    vh = Arnoldi{Float64}(8, 5)
-    V, H = vh.V, vh.H
-    V[:,1] .= 0.0
-    V[1,1] = 1.0
-
-    # Then {v, Av, A²v, A³v}
-    # is an invariant subspace
-    iterate_arnoldi!(A, vh, 1:5)
-
-    @test norm(V' * V - I) < 1e-10
-    @test iszero(H[5, 4])
+    for T in (Float64, BigFloat)
+        A = [rand(T,4,4)  zeros(T,4,4); 
+            zeros(T,4,4) rand(T,4,4)]
+
+        # and an initial vector [1; 0; ... 0]
+        vh = Arnoldi{T}(8, 5)
+        V, H = vh.V, vh.H
+        V[:,1] .= zero(T)
+        V[1,1] = one(T)
+
+        # Then {v, Av, A²v, A³v}
+        # is an invariant subspace
+        iterate_arnoldi!(A, vh, 1:5)
+
+        @test norm(V' * V - I) < sqrt(eps(T))/100
+        @test iszero(H[5, 4])
+    end
 end

test/givens_rotation.jl

@@ -6,7 +6,7 @@ using ArnoldiMethod: Hessenberg, Rotation2, Rotation3
 
 @testset "Givens rotation" begin
     @testset "Single rotation" begin
-        @testset "lmul!" for T in (Float64, ComplexF64)
+        @testset "lmul!" for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
             A = rand(T, 6, 5)
 
             G = Rotation2(rand(real(T)), rand(T), 2)
@@ -16,7 +16,7 @@ using ArnoldiMethod: Hessenberg, Rotation2, Rotation3
             @test G_mat * A ≈ lmul!(G, copy(A))
         end
 
-        @testset "rmul!" for T in (Float64, ComplexF64)
+        @testset "rmul!" for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
             A = rand(T, 10, 5)
 
             G = Rotation2(rand(real(T)), rand(T), 2)
@@ -28,7 +28,7 @@ using ArnoldiMethod: Hessenberg, Rotation2, Rotation3
     end
 
     @testset "Double rotation" begin
-        @testset "lmul!" for T in (Float64, ComplexF64)
+        @testset "lmul!" for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
             A = rand(T, 6, 5)
 
             G = Rotation3(rand(real(T)), rand(T), rand(real(T)), rand(T), 2)
@@ -38,7 +38,7 @@ using ArnoldiMethod: Hessenberg, Rotation2, Rotation3
             @test G_mat * A ≈ lmul!(G, copy(A))
         end
 
-        @testset "rmul!" for T in (Float64, ComplexF64)
+        @testset "rmul!" for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
             A = rand(T, 10, 5)
 
             G = Rotation3(rand(real(T)), rand(T), rand(real(T)), rand(T), 2)

test/partial_schur.jl

@@ -10,19 +10,20 @@ using LinearAlgebra
     # the stopping criterion works for these zero eigenvalues, we don't take
     # into account that it finds the complete eigenspace of the 0 eigenvalue.
     # (we have too few basis vectors for it anyways)
+    for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
+        A = rand(T,10, 3)
 
-    A = rand(10, 3)
+        # Rank 3 matrix.
+        B = A * A' 
 
-    # Rank 3 matrix.
-    B = A * A' 
+        schur, history = partialschur(B, nev = 5, mindim = 5, maxdim = 7, tol = eps())
 
-    schur, history = partialschur(B, nev = 5, mindim = 5, maxdim = 7, tol = eps())
-
-    @test history.converged
-    @test history.mvproducts == 7
-    @test norm(schur.Q'schur.Q - I) < 100eps()
-    @test norm(B * schur.Q - schur.Q * schur.R) < 100eps()
-    @test norm(diag(schur.R)[4:5]) < 100eps()
+        @test history.converged
+        @test history.mvproducts == 7
+        @test norm(schur.Q'schur.Q - I) < 100eps(real(T))
+        @test norm(B * schur.Q - schur.Q * schur.R) < 100eps(real(T))
+        @test norm(diag(schur.R)[4:5]) < 100eps(real(T))
+    end
 end
 
 @testset "Right number type" begin

test/runtests.jl

@@ -2,6 +2,7 @@ using Test
 using LinearAlgebra
 using SparseArrays
 using Random
+using GenericSchur
 
 include("expansion.jl")
 include("givens_rotation.jl")

test/schur_to_eigen.jl

@@ -2,11 +2,12 @@ using Test
 using LinearAlgebra, SparseArrays
 using ArnoldiMethod: partialschur, partialeigen
 using Random
+using GenericSchur
 
-@testset "Schur to eigen $T take $i" for T in (Float64,ComplexF64), i in 1:10
+@testset "Schur to eigen $T take $i" for T in (Float64,ComplexF64,BigFloat,Complex{BigFloat}), i in 1:10
     Random.seed!(i)
-    A = spdiagm(0 => 1:100) + sprand(100, 100, 0.01)
-    ε = 1e-7
+    A = spdiagm(0 => 1:100) + sprand(T,100, 100, 0.01)
+    ε = sqrt(eps(real(T)))
     minim, maxim = 10, 20
 
     decomp, history = partialschur(A, nev=minim, tol=ε, restarts=200)
@@ -18,4 +19,4 @@ using Random
     for i = 1 : minim
         @test norm(A * vecs[:,i] - vecs[:,i] * vals[i]) < ε * abs(vals[i])
     end
-end
+end