From 7a27157cef88c4c5c1f747f0ce5f53cb1ce08f48 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 7 Mar 2024 18:21:19 +0100 Subject: [PATCH] Refactor tests II --- test/eigsolve.jl | 285 +++++++++++++----------------------------- test/expintegrator.jl | 196 ++++++++--------------------- test/factorize.jl | 194 ++++++++-------------------- test/geneigsolve.jl | 103 +++++---------- test/gklfactorize.jl | 179 +++++--------------------- test/linsolve.jl | 193 ++++++++-------------------- test/schursolve.jl | 164 +++++++----------------- test/svdsolve.jl | 114 ++++------------- test/testsetup.jl | 43 ++++++- 9 files changed, 413 insertions(+), 1058 deletions(-) diff --git a/test/eigsolve.jl b/test/eigsolve.jl index e55f96c..717d95e 100644 --- a/test/eigsolve.jl +++ b/test/eigsolve.jl @@ -1,72 +1,57 @@ -@testset "Lanczos - eigsolve full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Lanczos - eigsolve full ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) .- one(T) / 2 A = (A + A') / 2 v = rand(T, (n,)) n1 = div(n, 2) - D1, V1, info = eigsolve(A, v, n1, :SR; orth=orth, krylovdim=n, + D1, V1, info = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n1, :SR; + krylovdim=n, maxiter=1, tol=tolerance(T), verbosity=2) - @test KrylovKit.eigselector(A, scalartype(v); orth=orth, krylovdim=n, maxiter=1, - tol=tolerance(T)) isa Lanczos + @test KrylovKit.eigselector(wrapop(A, Val(mode)), scalartype(v); krylovdim=n, + maxiter=1, + tol=tolerance(T), ishermitian=true) isa Lanczos n2 = n - n1 - alg = Lanczos(; orth=orth, krylovdim=2 * n, maxiter=1, tol=tolerance(T), + alg = Lanczos(; krylovdim=2 * n, maxiter=1, tol=tolerance(T), verbosity=1) - D2, V2, info = @constinferred eigsolve(A, v, n2, :LR, alg) + D2, V2, info = @constinferred eigsolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), + n2, :LR, alg) @test vcat(D1[1:n1], reverse(D2[1:n2])) ≊ eigvals(A) - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1' * U1 ≈ I @test U2' * U2 ≈ I @test A * U1 ≈ U1 * Diagonal(D1) @test A * U2 ≈ U2 * Diagonal(D2) - _ = eigsolve(A, v, n + 1, :LM; orth=orth, krylovdim=2n, + _ = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n + 1, :LM; + krylovdim=2n, maxiter=1, tol=tolerance(T), verbosity=0) end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - A = (A + A') / 2 - v = MinimalVec{IP}(rand(T, (n,))) - n1 = div(n, 2) - D1, V1, info = eigsolve(wrapop(A), v, n1, :SR; krylovdim=n, - maxiter=1, tol=tolerance(T), verbosity=2) - @test KrylovKit.eigselector(wrapop(A), scalartype(v); krylovdim=n, maxiter=1, - tol=tolerance(T), ishermitian=true) isa Lanczos - n2 = n - n1 - alg = Lanczos(; krylovdim=2 * n, maxiter=1, tol=tolerance(T), - verbosity=1) - D2, V2, info = @constinferred eigsolve(wrapop(A), v, n2, :LR, alg) - @test vcat(D1[1:n1], reverse(D2[1:n2])) ≊ eigvals(A) - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1' * U1 ≈ I - @test U2' * U2 ≈ I - - @test A * U1 ≈ U1 * Diagonal(D1) - @test A * U2 ≈ U2 * Diagonal(D2) - - _ = eigsolve(wrapop(A), v, n + 1, :LM; krylovdim=2n, - maxiter=1, tol=tolerance(T), verbosity=0) - end end -@testset "Lanczos - eigsolve iteratively" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Lanczos - eigsolve iteratively ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (N, N)) .- one(T) / 2 A = (A + A') / 2 v = rand(T, (N,)) - alg = Lanczos(; orth=orth, krylovdim=2 * n, maxiter=10, + alg = Lanczos(; krylovdim=2 * n, maxiter=10, tol=tolerance(T), eager=true) - D1, V1, info1 = @constinferred eigsolve(A, v, n, :SR, alg) - D2, V2, info2 = eigsolve(A, v, n, :LR, alg) + D1, V1, info1 = @constinferred eigsolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), n, :SR, alg) + D2, V2, info2 = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, :LR, + alg) l1 = info1.converged l2 = info2.converged @@ -75,140 +60,90 @@ end @test D1[1:l1] ≈ eigvals(A)[1:l1] @test D2[1:l2] ≈ eigvals(A)[N:-1:(N - l2 + 1)] - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1' * U1 ≈ I @test U2' * U2 ≈ I - R1 = stack(info1.residual) - R2 = stack(info2.residual) + R1 = stack(unwrapvec, info1.residual) + R2 = stack(unwrapvec, info2.residual) @test A * U1 ≈ U1 * Diagonal(D1) + R1 @test A * U2 ≈ U2 * Diagonal(D2) + R2 end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - A = (A + A') / 2 - v = MinimalVec{IP}(rand(T, (N,))) - alg = Lanczos(; krylovdim=2 * n, maxiter=10, - tol=tolerance(T), eager=true) - D1, V1, info1 = @constinferred eigsolve(wrapop(A), v, n, :SR, alg) - D2, V2, info2 = eigsolve(wrapop(A), v, n, :LR, alg) - - l1 = info1.converged - l2 = info2.converged - @test l1 > 0 - @test l2 > 0 - @test D1[1:l1] ≈ eigvals(A)[1:l1] - @test D2[1:l2] ≈ eigvals(A)[N:-1:(N - l2 + 1)] - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1' * U1 ≈ I - @test U2' * U2 ≈ I - - R1 = stack(unwrap, info1.residual) - R2 = stack(unwrap, info2.residual) - @test A * U1 ≈ U1 * Diagonal(D1) + R1 - @test A * U2 ≈ U2 * Diagonal(D2) + R2 - end end -@testset "Arnoldi - eigsolve full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Arnoldi - eigsolve full ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) .- one(T) / 2 v = rand(T, (n,)) n1 = div(n, 2) - D1, V1, info1 = eigsolve(A, v, n1, :SR; orth=orth, krylovdim=n, + D1, V1, info1 = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n1, :SR; + orth=orth, krylovdim=n, maxiter=1, tol=tolerance(T), verbosity=2) - @test KrylovKit.eigselector(A, eltype(v); orth=orth, krylovdim=n, maxiter=1, + @test KrylovKit.eigselector(wrapop(A, Val(mode)), eltype(v); orth=orth, + krylovdim=n, maxiter=1, tol=tolerance(T)) isa Arnoldi n2 = n - n1 alg = Arnoldi(; orth=orth, krylovdim=2 * n, maxiter=1, tol=tolerance(T), verbosity=1) - D2, V2, info2 = @constinferred eigsolve(A, v, n2, :LR, alg) + D2, V2, info2 = @constinferred eigsolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), n2, :LR, alg) D = sort(sort(eigvals(A); by=imag, rev=true); alg=MergeSort, by=real) D2′ = sort(sort(D2; by=imag, rev=true); alg=MergeSort, by=real) @test vcat(D1[1:n1], D2′[(end - n2 + 1):end]) ≈ D - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test A * U1 ≈ U1 * Diagonal(D1) @test A * U2 ≈ U2 * Diagonal(D2) if T <: Complex n1 = div(n, 2) - D1, V1, info = eigsolve(A, v, n1, :SI, alg) + D1, V1, info = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n1, + :SI, + alg) n2 = n - n1 - D2, V2, info = eigsolve(A, v, n2, :LI, alg) + D2, V2, info = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n2, + :LI, + alg) D = sort(eigvals(A); by=imag) @test vcat(D1[1:n1], reverse(D2[1:n2])) ≊ D - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test A * U1 ≈ U1 * Diagonal(D1) @test A * U2 ≈ U2 * Diagonal(D2) end - _ = eigsolve(A, v, n + 1, :LM; orth=orth, krylovdim=2n, + _ = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n + 1, :LM; orth=orth, + krylovdim=2n, maxiter=1, tol=tolerance(T), verbosity=0) end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - v = MinimalVec{IP}(rand(T, (n,))) - n1 = div(n, 2) - D1, V1, info1 = eigsolve(wrapop(A), v, n1, :SR; krylovdim=n, - maxiter=1, tol=tolerance(T), verbosity=2) - @test KrylovKit.eigselector(wrapop(A), scalartype(v); krylovdim=n, maxiter=1, - tol=tolerance(T), ishermitian=false) isa Arnoldi - n2 = n - n1 - alg = Arnoldi(; krylovdim=2 * n, maxiter=1, tol=tolerance(T), - verbosity=1) - D2, V2, info2 = @constinferred eigsolve(wrapop(A), v, n2, :LR, alg) - D = sort(sort(eigvals(A); by=imag, rev=true); alg=MergeSort, by=real) - D2′ = sort(sort(D2; by=imag, rev=true); alg=MergeSort, by=real) - @test vcat(D1[1:n1], D2′[(end - n2 + 1):end]) ≈ D - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test A * U1 ≈ U1 * Diagonal(D1) - @test A * U2 ≈ U2 * Diagonal(D2) - - n1 = div(n, 2) - D1, V1, info = eigsolve(wrapop(A), v, n1, :SI, alg) - n2 = n - n1 - D2, V2, info = eigsolve(wrapop(A), v, n2, :LI, alg) - D = sort(eigvals(A); by=imag) - - @test vcat(D1[1:n1], reverse(D2[1:n2])) ≊ D - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test A * U1 ≈ U1 * Diagonal(D1) - @test A * U2 ≈ U2 * Diagonal(D2) - - _ = eigsolve(wrapop(A), v, n + 1, :LM; krylovdim=2n, - maxiter=1, tol=tolerance(T), verbosity=0) - end end -@testset "Arnoldi - eigsolve iteratively" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Arnoldi - eigsolve iteratively ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (N, N)) .- one(T) / 2 v = rand(T, (N,)) - alg = Arnoldi(; orth=orth, krylovdim=3 * n, maxiter=20, + alg = Arnoldi(; krylovdim=3 * n, maxiter=20, tol=tolerance(T), eager=true) - D1, V1, info1 = @constinferred eigsolve(A, v, n, :SR, alg) - D2, V2, info2 = eigsolve(A, v, n, :LR, alg) - D3, V3, info3 = eigsolve(A, v, n, :LM, alg) + D1, V1, info1 = @constinferred eigsolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), n, :SR, alg) + D2, V2, info2 = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, :LR, + alg) + D3, V3, info3 = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, :LM, + alg) D = sort(eigvals(A); by=imag, rev=true) l1 = info1.converged @@ -223,19 +158,21 @@ end # in absolute value, so we perform a second sort afterwards using the real part @test D3[1:l3] ≊ sort(D; by=abs, rev=true)[1:l3] - U1 = stack(V1) - U2 = stack(V2) - U3 = stack(V3) - R1 = stack(info1.residual) - R2 = stack(info2.residual) - R3 = stack(info3.residual) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) + U3 = stack(unwrapvec, V3) + R1 = stack(unwrapvec, info1.residual) + R2 = stack(unwrapvec, info2.residual) + R3 = stack(unwrapvec, info3.residual) @test A * U1 ≈ U1 * Diagonal(D1) + R1 @test A * U2 ≈ U2 * Diagonal(D2) + R2 @test A * U3 ≈ U3 * Diagonal(D3) + R3 if T <: Complex - D1, V1, info1 = eigsolve(A, v, n, :SI, alg) - D2, V2, info2 = eigsolve(A, v, n, :LI, alg) + D1, V1, info1 = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, + :SI, alg) + D2, V2, info2 = eigsolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, + :LI, alg) D = eigvals(A) l1 = info1.converged @@ -245,65 +182,13 @@ end @test D1[1:l1] ≈ sort(D; by=imag)[1:l1] @test D2[1:l2] ≈ sort(D; by=imag, rev=true)[1:l2] - U1 = stack(V1) - U2 = stack(V2) - R1 = stack(info1.residual) - R2 = stack(info2.residual) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) + R1 = stack(unwrapvec, info1.residual) + R2 = stack(unwrapvec, info2.residual) @test A * U1 ≈ U1 * Diagonal(D1) + R1 @test A * U2 ≈ U2 * Diagonal(D2) + R2 end end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - v = MinimalVec{IP}(rand(T, (N,))) - alg = Arnoldi(; krylovdim=3 * n, maxiter=20, - tol=tolerance(T), eager=true) - D1, V1, info1 = @constinferred eigsolve(wrapop(A), v, n, :SR, alg) - D2, V2, info2 = eigsolve(wrapop(A), v, n, :LR, alg) - D3, V3, info3 = eigsolve(wrapop(A), v, n, :LM, alg) - D = sort(eigvals(A); by=imag, rev=true) - - l1 = info1.converged - l2 = info2.converged - l3 = info3.converged - @test l1 > 0 - @test l2 > 0 - @test l3 > 0 - @test D1[1:l1] ≊ sort(D; alg=MergeSort, by=real)[1:l1] - @test D2[1:l2] ≊ sort(D; alg=MergeSort, by=real, rev=true)[1:l2] - # sorting by abs does not seem very reliable if two distinct eigenvalues are close - # in absolute value, so we perform a second sort afterwards using the real part - @test D3[1:l3] ≊ sort(D; by=abs, rev=true)[1:l3] - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - U3 = stack(unwrap, V3) - R1 = stack(unwrap, info1.residual) - R2 = stack(unwrap, info2.residual) - R3 = stack(unwrap, info3.residual) - @test A * U1 ≈ U1 * Diagonal(D1) + R1 - @test A * U2 ≈ U2 * Diagonal(D2) + R2 - @test A * U3 ≈ U3 * Diagonal(D3) + R3 - - D1, V1, info1 = eigsolve(wrapop(A), v, n, :SI, alg) - D2, V2, info2 = eigsolve(wrapop(A), v, n, :LI, alg) - D = eigvals(A) - - l1 = info1.converged - l2 = info2.converged - @test l1 > 0 - @test l2 > 0 - @test D1[1:l1] ≈ sort(D; by=imag)[1:l1] - @test D2[1:l2] ≈ sort(D; by=imag, rev=true)[1:l2] - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - R1 = stack(unwrap, info1.residual) - R2 = stack(unwrap, info2.residual) - @test A * U1 ≈ U1 * Diagonal(D1) + R1 - @test A * U2 ≈ U2 * Diagonal(D2) + R2 - end end diff --git a/test/expintegrator.jl b/test/expintegrator.jl index 9946321..fd48972 100644 --- a/test/expintegrator.jl +++ b/test/expintegrator.jl @@ -12,9 +12,12 @@ function ϕ(A, v, p) return exp(A′)[1:m, end] end -@testset "Lanczos - expintegrator full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Lanczos - expintegrator full ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) .- one(T) / 2 A = (A + A') / 2 V = one(A) @@ -22,7 +25,11 @@ end alg = Lanczos(; orth=orth, krylovdim=n, maxiter=2, tol=tolerance(T), verbosity=2) for k in 1:n - W[:, k] = first(@constinferred exponentiate(A, 1, view(V, :, k), alg)) + W[:, k] = unwrapvec(first(@constinferred exponentiate(wrapop(A, Val(mode)), + 1, + wrapvec(view(V, :, k), + Val(mode)), + alg))) end @test W ≈ exp(A) @@ -33,60 +40,37 @@ end randn(real(T)) + im * randn(real(T))) for p in 1:pmax u = ntuple(i -> rand(T, n), p + 1) - w, info = @constinferred expintegrator(A, t, u, alg) + w, info = @constinferred expintegrator(wrapop(A, Val(mode)), t, + wrapvec.(u, Ref(Val(mode))), alg) w2 = exp(t * A) * u[1] for j in 1:p w2 .+= t^j * ϕ(t * A, u[j + 1], j) end @test info.converged > 0 - @test w2 ≈ w + @test w2 ≈ unwrapvec(w) end end end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - A = (A + A') / 2 - V = one(A) - W = zero(A) - alg = Lanczos(; krylovdim=n, maxiter=2, tol=tolerance(T), verbosity=2) - for k in 1:n - W[:, k] = unwrap(first(@constinferred exponentiate(wrapop(A), 1, - MinimalVec{IP}(view(V, :, k)), - alg))) - end - @test W ≈ exp(A) - - pmax = 5 - alg = Lanczos(; krylovdim=n, maxiter=2, tol=tolerance(T), verbosity=1) - for t in (rand(real(T)), -rand(real(T)), im * randn(real(T)), - randn(real(T)) + im * randn(real(T))) - for p in 1:pmax - u = ntuple(i -> MinimalVec{IP}(rand(T, n)), p + 1) - w, info = @constinferred expintegrator(wrapop(A), t, u, alg) - w2 = exp(t * A) * unwrap(u[1]) - for j in 1:p - w2 .+= t^j * ϕ(t * A, unwrap(u[j + 1]), j) - end - @test info.converged > 0 - @test w2 ≈ unwrap(w) - end - end - end end -@testset "Arnoldi - expintegrator full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Arnoldi - expintegrator full ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) .- one(T) / 2 V = one(A) W = zero(A) alg = Arnoldi(; orth=orth, krylovdim=n, maxiter=2, tol=tolerance(T), verbosity=2) for k in 1:n - W[:, k] = first(@constinferred exponentiate(A, 1, view(V, :, k), alg)) + W[:, k] = unwrapvec(first(@constinferred exponentiate(wrapop(A, Val(mode)), + 1, + wrapvec(view(V, :, k), + Val(mode)), + alg))) end @test W ≈ exp(A) @@ -97,52 +81,27 @@ end randn(real(T)) + im * randn(real(T))) for p in 1:pmax u = ntuple(i -> rand(T, n), p + 1) - w, info = @constinferred expintegrator(A, t, u, alg) + w, info = @constinferred expintegrator(wrapop(A, Val(mode)), t, + wrapvec.(u, Ref(Val(mode))), alg) w2 = exp(t * A) * u[1] for j in 1:p w2 .+= t^j * ϕ(t * A, u[j + 1], j) end @test info.converged > 0 - @test w2 ≈ w - end - end - end - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - V = one(A) - W = zero(A) - alg = Arnoldi(; krylovdim=n, maxiter=2, tol=tolerance(T), verbosity=2) - for k in 1:n - W[:, k] = unwrap(first(@constinferred exponentiate(wrapop(A), 1, - MinimalVec{IP}(view(V, :, k)), - alg))) - end - @test W ≈ exp(A) - - pmax = 5 - alg = Arnoldi(; krylovdim=n, maxiter=2, tol=tolerance(T), verbosity=1) - for t in (rand(real(T)), -rand(real(T)), im * randn(real(T)), - randn(real(T)) + im * randn(real(T))) - for p in 1:pmax - u = ntuple(i -> MinimalVec{IP}(rand(T, n)), p + 1) - w, info = @constinferred expintegrator(wrapop(A), t, u, alg) - w2 = exp(t * A) * unwrap(u[1]) - for j in 1:p - w2 .+= t^j * ϕ(t * A, unwrap(u[j + 1]), j) + @test w2 ≈ unwrapvec(w) end - @test info.converged > 0 - @test w2 ≈ unwrap(w) end end end end -@testset "Lanczos - expintegrator iteratively" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Lanczos - expintegrator iteratively ($mode)" for mode in + (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (N, N)) .- one(T) / 2 A = (A + A') / 2 s = norm(eigvals(A), 1) @@ -152,7 +111,8 @@ end randn(real(T)) + im * randn(real(T))) for p in 1:pmax u = ntuple(i -> rand(T, N), p + 1) - w1, info = @constinferred expintegrator(A, t, u...; + w1, info = @constinferred expintegrator(wrapop(A, Val(mode)), t, + wrapvec.(u, Ref(Val(mode)))...; maxiter=100, krylovdim=n, eager=true) @assert info.converged > 0 @@ -160,48 +120,25 @@ end for j in 1:p w2 .+= t^j * ϕ(t * A, u[j + 1], j) end - @test w2 ≈ w1 - w1, info = @constinferred expintegrator(A, t, u...; + @test w2 ≈ unwrapvec(w1) + w1, info = @constinferred expintegrator(wrapop(A, Val(mode)), t, + wrapvec.(u, Ref(Val(mode)))...; maxiter=100, krylovdim=n, tol=1e-3, eager=true) - @test w1 ≈ w2 atol = 1e-2 * abs(t) - end - end - end - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - A = (A + A') / 2 - s = norm(eigvals(A), 1) - rmul!(A, 1 / (10 * s)) - pmax = 5 - for t in (rand(real(T)), -rand(real(T)), im * randn(real(T)), - randn(real(T)) + im * randn(real(T))) - for p in 1:pmax - u = ntuple(i -> MinimalVec{IP}(rand(T, N)), p + 1) - w1, info = @constinferred expintegrator(wrapop(A), t, u...; - maxiter=100, krylovdim=n, - eager=true) - @assert info.converged > 0 - w2 = exp(t * A) * unwrap(u[1]) - for j in 1:p - w2 .+= t^j * ϕ(t * A, unwrap(u[j + 1]), j) + @test unwrapvec(w1) ≈ w2 atol = 1e-2 * abs(t) end - @test w2 ≈ unwrap(w1) - w1, info = @constinferred expintegrator(wrapop(A), t, u...; - maxiter=100, krylovdim=n, - tol=1e-3, eager=true) - @test unwrap(w1) ≈ w2 atol = 1e-2 * abs(t) end end end end -@testset "Arnoldi - expintegrator iteratively" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Arnoldi - expintegrator iteratively ($mode)" for mode in + (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (N, N)) .- one(T) / 2 s = norm(eigvals(A), 1) rmul!(A, 1 / (10 * s)) @@ -210,7 +147,8 @@ end randn(real(T)) + im * randn(real(T))) for p in 1:pmax u = ntuple(i -> rand(T, N), p + 1) - w1, info = @constinferred expintegrator(A, t, u...; + w1, info = @constinferred expintegrator(wrapop(A, Val(mode)), t, + wrapvec.(u, Ref(Val(mode)))...; maxiter=100, krylovdim=n, eager=true) @test info.converged > 0 @@ -218,39 +156,13 @@ end for j in 1:p w2 .+= t^j * ϕ(t * A, u[j + 1], j) end - @test w2 ≈ w1 - w1, info = @constinferred expintegrator(A, t, u...; + @test w2 ≈ unwrapvec(w1) + w1, info = @constinferred expintegrator(wrapop(A, Val(mode)), t, + wrapvec.(u, Ref(Val(mode)))...; maxiter=100, krylovdim=n, tol=1e-3, eager=true) - @test w1 ≈ w2 atol = 1e-2 * abs(t) - end - end - end - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - s = norm(eigvals(A), 1) - rmul!(A, 1 / (10 * s)) - pmax = 5 - for t in (rand(real(T)), -rand(real(T)), im * randn(real(T)), - randn(real(T)) + im * randn(real(T))) - for p in 1:pmax - u = ntuple(i -> MinimalVec{IP}(rand(T, N)), p + 1) - w1, info = @constinferred expintegrator(wrapop(A), t, u...; - maxiter=100, krylovdim=n, - eager=true) - @test info.converged > 0 - w2 = exp(t * A) * unwrap(u[1]) - for j in 1:p - w2 .+= t^j * ϕ(t * A, unwrap(u[j + 1]), j) + @test unwrapvec(w1) ≈ w2 atol = 1e-2 * abs(t) end - @test w2 ≈ unwrap(w1) - w1, info = @constinferred expintegrator(wrapop(A), t, u...; - maxiter=100, krylovdim=n, - tol=1e-3, eager=true) - @test unwrap(w1) ≈ w2 atol = 1e-2 * abs(t) end end end diff --git a/test/factorize.jl b/test/factorize.jl index 2069788..467b195 100644 --- a/test/factorize.jl +++ b/test/factorize.jl @@ -1,15 +1,15 @@ # Test complete Lanczos factorization -wrapop(A) = function (v::MinimalVec) - return MinimalVec{isinplace(v)}(A * unwrap(v)) -end +@testset "Complete Lanczos factorization ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (cgs2,) -@testset "Complete Lanczos factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) # tests fail miserably for cgs and mgs + @testset for T in scalartypes + @testset for orth in orths # tests fail miserably for cgs and mgs A = rand(T, (n, n)) v = rand(T, (n,)) A = (A + A') - iter = LanczosIterator(A, v, orth) + iter = LanczosIterator(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), orth) verbosity = 1 fact = @constinferred initialize(iter; verbosity=verbosity) while length(fact) < n @@ -17,7 +17,7 @@ end verbosity = 0 end - V = stack(basis(fact)) + V = stack(unwrapvec, basis(fact)) H = rayleighquotient(fact) @test normres(fact) < 10 * n * eps(real(T)) @test V' * V ≈ I @@ -28,40 +28,19 @@ end @test rayleighquotient(last(states)) ≈ H end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) - A += A' - - v = MinimalVec{IP}(rand(T, (n,))) - iter = LanczosIterator(wrapop(A), v) - verbosity = 1 - fact = @constinferred initialize(iter; verbosity=verbosity) - while length(fact) < n - @constinferred expand!(iter, fact; verbosity=verbosity) - verbosity = 0 - end - - V = stack(unwrap, basis(fact)) - H = rayleighquotient(fact) - @test normres(fact) < 10 * n * eps(real(T)) - @test V' * V ≈ I - @test A * V ≈ V * H - - @constinferred initialize!(iter, deepcopy(fact); verbosity=1) - states = collect(Iterators.take(iter, n)) # collect tests size and eltype? - @test rayleighquotient(last(states)) ≈ H - end end # Test complete Arnoldi factorization -@testset "Complete Arnoldi factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs, mgs, cgs2, mgs2, cgsr, mgsr) +@testset "Complete Arnoldi factorization ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs, mgs, cgs2, mgs2, cgsr, mgsr) : (cgs2,) + + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) v = rand(T, (n,)) - iter = ArnoldiIterator(A, v, orth) + iter = ArnoldiIterator(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), orth) verbosity = 1 fact = @constinferred initialize(iter; verbosity=verbosity) while length(fact) < n @@ -69,7 +48,7 @@ end verbosity = 0 end - V = stack(basis(fact)) + V = stack(unwrapvec, basis(fact)) H = rayleighquotient(fact) factor = (orth == cgs || orth == mgs ? 250 : 10) @test normres(fact) < factor * n * eps(real(T)) @@ -81,37 +60,18 @@ end @test rayleighquotient(last(states)) ≈ H end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) - - v = MinimalVec{IP}(rand(T, (n,))) - iter = ArnoldiIterator(wrapop(A), v) - verbosity = 1 - fact = @constinferred initialize(iter; verbosity=verbosity) - while length(fact) < n - @constinferred expand!(iter, fact; verbosity=verbosity) - verbosity = 0 - end - - V = stack(unwrap, basis(fact)) - H = rayleighquotient(fact) - @test normres(fact) < 10 * n * eps(real(T)) - @test V' * V ≈ I - @test A * V ≈ V * H - - @constinferred initialize!(iter, deepcopy(fact); verbosity=1) - states = collect(Iterators.take(iter, n)) # collect tests size and eltype? - @test rayleighquotient(last(states)) ≈ H - end end # Test incomplete Lanczos factorization -@testset "Incomplete Lanczos factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int}) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) # tests fail miserably for cgs and mgs - if T == Complex{Int} +@testset "Incomplete Lanczos factorization ($mode)" for mode in + (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? + (Float32, Float64, ComplexF32, ComplexF64, Complex{Int}) : (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (cgs2,) + + @testset for T in scalartypes + @testset for orth in orths # tests fail miserably for cgs and mgs + if T === Complex{Int} A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N)) v = rand(-100:100, (N,)) else @@ -119,24 +79,25 @@ end v = rand(T, (N,)) end A = (A + A') - iter = @constinferred LanczosIterator(A, v, orth) + iter = @constinferred LanczosIterator(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), + orth) krylovdim = n fact = @constinferred initialize(iter) while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim @constinferred expand!(iter, fact) - - Ṽ, H, r, β, e = fact - V = stack(Ṽ) + Ṽ, H, r̃, β, e = fact + V = stack(unwrapvec, Ṽ) + r = unwrapvec(r̃) @test V' * V ≈ I @test norm(r) ≈ β @test A * V ≈ V * H + r * e' end fact = @constinferred shrink!(fact, div(n, 2)) - B = @constinferred basis(fact) - V = stack(B) + V = stack(unwrapvec, @constinferred basis(fact)) H = @constinferred rayleighquotient(fact) - r = @constinferred residual(fact) + r = @constinferred unwrapvec(residual(fact)) β = @constinferred normres(fact) e = @constinferred rayleighextension(fact) @test V' * V ≈ I @@ -144,67 +105,42 @@ end @test A * V ≈ V * H + r * e' end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) - A += A' - - v = MinimalVec{IP}(rand(T, (N,))) - iter = @constinferred LanczosIterator(wrapop(A), v) - krylovdim = n - fact = @constinferred initialize(iter) - while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim - @constinferred expand!(iter, fact) - - Ṽ, H, r̃, β, e = fact - V = stack(unwrap, Ṽ) - r = unwrap(r̃) - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ V * H + r * e' - end - - fact = @constinferred shrink!(fact, div(n, 2)) - V = stack(unwrap, @constinferred basis(fact)) - H = @constinferred rayleighquotient(fact) - r = unwrap(@constinferred residual(fact)) - β = @constinferred normres(fact) - e = @constinferred rayleighextension(fact) - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ V * H + r * e' - end end # Test incomplete Arnoldi factorization -@testset "Incomplete Arnoldi factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int}) - @testset for orth in (cgs, mgs, cgs2, mgs2, cgsr, mgsr) - if T == Complex{Int} +@testset "Incomplete Arnoldi factorization ($mode)" for mode in + (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? + (Float32, Float64, ComplexF32, ComplexF64, Complex{Int}) : (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (cgs2,) + + @testset for T in scalartypes + @testset for orth in orths + if T === Complex{Int} A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N)) v = rand(-100:100, (N,)) else A = rand(T, (N, N)) v = rand(T, (N,)) end - iter = @constinferred ArnoldiIterator(A, v, orth) + iter = @constinferred ArnoldiIterator(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), orth) krylovdim = 3 * n fact = @constinferred initialize(iter) while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim @constinferred expand!(iter, fact) - - Ṽ, H, r, β, e = fact - V = stack(Ṽ) + Ṽ, H, r̃, β, e = fact + V = stack(unwrapvec, Ṽ) + r = unwrapvec(r̃) @test V' * V ≈ I @test norm(r) ≈ β @test A * V ≈ V * H + r * e' end fact = @constinferred shrink!(fact, div(n, 2)) - V = stack(@constinferred basis(fact)) + V = stack(unwrapvec, @constinferred basis(fact)) H = @constinferred rayleighquotient(fact) - r = @constinferred residual(fact) + r = unwrapvec(@constinferred residual(fact)) β = @constinferred normres(fact) e = @constinferred rayleighextension(fact) @test V' * V ≈ I @@ -212,34 +148,4 @@ end @test A * V ≈ V * H + r * e' end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) - - v = MinimalVec{IP}(rand(T, (N,))) - iter = @constinferred ArnoldiIterator(wrapop(A), v) - krylovdim = 3 * n - fact = @constinferred initialize(iter) - while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim - @constinferred expand!(iter, fact) - - Ṽ, H, r̃, β, e = fact - V = stack(unwrap, Ṽ) - r = unwrap(r̃) - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ V * H + r * e' - end - - fact = @constinferred shrink!(fact, div(n, 2)) - V = stack(unwrap, @constinferred basis(fact)) - H = @constinferred rayleighquotient(fact) - r = unwrap(@constinferred residual(fact)) - β = @constinferred normres(fact) - e = @constinferred rayleighextension(fact) - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ V * H + r * e' - end end diff --git a/test/geneigsolve.jl b/test/geneigsolve.jl index 9105314..d154b81 100644 --- a/test/geneigsolve.jl +++ b/test/geneigsolve.jl @@ -1,7 +1,10 @@ -@testset "GolubYe - geneigsolve full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "GolubYe - geneigsolve full ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) .- one(T) / 2 A = (A + A') / 2 B = rand(T, (n, n)) .- one(T) / 2 @@ -10,21 +13,26 @@ alg = GolubYe(; orth=orth, krylovdim=n, maxiter=1, tol=tolerance(T), verbosity=1) n1 = div(n, 2) - D1, V1, info = @constinferred geneigsolve((A, B), v, + D1, V1, info = @constinferred geneigsolve((wrapop(A, Val(mode)), + wrapop(B, Val(mode))), + wrapvec(v, Val(mode)), n1, :SR; orth=orth, krylovdim=n, maxiter=1, tol=tolerance(T), ishermitian=true, isposdef=true, verbosity=2) - @test KrylovKit.geneigselector((A, B), scalartype(v); orth=orth, krylovdim=n, + @test KrylovKit.geneigselector((wrapop(A, Val(mode)), wrapop(B, Val(mode))), + scalartype(v); orth=orth, krylovdim=n, maxiter=1, tol=tolerance(T), ishermitian=true, isposdef=true) isa GolubYe n2 = n - n1 - D2, V2, info = @constinferred geneigsolve((A, B), v, + D2, V2, info = @constinferred geneigsolve((wrapop(A, Val(mode)), + wrapop(B, Val(mode))), + wrapvec(v, Val(mode)), n2, :LR, alg) @test vcat(D1[1:n1], reverse(D2[1:n2])) ≈ eigvals(A, B) - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1' * B * U1 ≈ I @test U2' * B * U2 ≈ I @@ -32,35 +40,15 @@ @test A * U2 ≈ B * U2 * Diagonal(D2) end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - A = (A + A') / 2 - B = rand(T, (n, n)) .- one(T) / 2 - B = sqrt(B * B') - v = MinimalVec{IP}(rand(T, (n,))) - alg = GolubYe(; krylovdim=n, maxiter=1, tol=tolerance(T), verbosity=1) - n1 = div(n, 2) - D1, V1, info = @constinferred geneigsolve((wrapop(A), wrapop(B)), v, - n1, :SR, alg) - n2 = n - n1 - D2, V2, info = geneigsolve((wrapop(A), wrapop(B)), v, n2, :LR, alg) - @test vcat(D1[1:n1], reverse(D2[1:n2])) ≈ eigvals(A, B) - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1' * B * U1 ≈ I - @test U2' * B * U2 ≈ I - - @test A * U1 ≈ B * U1 * Diagonal(D1) - @test A * U2 ≈ B * U2 * Diagonal(D2) - end end -@testset "GolubYe - geneigsolve iteratively" begin - @testset for T in (Float64, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "GolubYe - geneigsolve iteratively ($mode)" for mode in + (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float64, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (N, N)) .- one(T) / 2 A = (A + A') / 2 B = rand(T, (N, N)) .- one(T) / 2 @@ -68,9 +56,12 @@ end v = rand(T, (N,)) alg = GolubYe(; orth=orth, krylovdim=3 * n, maxiter=100, tol=cond(B) * tolerance(T)) - D1, V1, info1 = @constinferred geneigsolve((A, B), v, + D1, V1, info1 = @constinferred geneigsolve((wrapop(A, Val(mode)), + wrapop(B, Val(mode))), + wrapvec(v, Val(mode)), n, :SR, alg) - D2, V2, info2 = geneigsolve((A, B), v, n, :LR, alg) + D2, V2, info2 = geneigsolve((wrapop(A, Val(mode)), wrapop(B, Val(mode))), + wrapvec(v, Val(mode)), n, :LR, alg) l1 = info1.converged l2 = info2.converged @@ -79,45 +70,15 @@ end @test D1[1:l1] ≊ eigvals(A, B)[1:l1] @test D2[1:l2] ≊ eigvals(A, B)[N:-1:(N - l2 + 1)] - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1' * B * U1 ≈ I @test U2' * B * U2 ≈ I - R1 = stack(info1.residual) - R2 = stack(info2.residual) + R1 = stack(unwrapvec, info1.residual) + R2 = stack(unwrapvec, info2.residual) @test A * U1 ≈ B * U1 * Diagonal(D1) + R1 @test A * U2 ≈ B * U2 * Diagonal(D2) + R2 end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - A = (A + A') / 2 - B = rand(T, (N, N)) .- one(T) / 2 - B = sqrt(B * B') - v = MinimalVec{IP}(rand(T, (N,))) - alg = GolubYe(; krylovdim=3 * n, maxiter=100, - tol=cond(B) * tolerance(T)) - D1, V1, info1 = geneigsolve((wrapop(A), wrapop(B)), v, n, :SR, alg) - D2, V2, info2 = geneigsolve((wrapop(A), wrapop(B)), v, n, :LR, alg) - - l1 = info1.converged - l2 = info2.converged - @test l1 > 0 - @test l2 > 0 - @test D1[1:l1] ≊ eigvals(A, B)[1:l1] - @test D2[1:l2] ≊ eigvals(A, B)[N:-1:(N - l2 + 1)] - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1' * B * U1 ≈ I - @test U2' * B * U2 ≈ I - - R1 = stack(unwrap, info1.residual) - R2 = stack(unwrap, info2.residual) - @test A * U1 ≈ B * U1 * Diagonal(D1) + R1 - @test A * U2 ≈ B * U2 * Diagonal(D2) + R2 - end end diff --git a/test/gklfactorize.jl b/test/gklfactorize.jl index 7d89a12..471e01b 100644 --- a/test/gklfactorize.jl +++ b/test/gklfactorize.jl @@ -1,23 +1,15 @@ -wrapop(A) = function (v, flag=Val(false)) - if flag === Val(true) - return MinimalVec{isinplace(v)}(A' * unwrap(v)) - else - return MinimalVec{isinplace(v)}(A * unwrap(v)) - end -end - -function wrapop2(A) - return (x -> MinimalVec{false}(A * unwrap(x)), - y -> MinimalVec{true}(A' * unwrap(y))) -end - # Test complete Golub-Kahan-Lanczos factorization -@testset "Complete Golub-Kahan-Lanczos factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Complete Golub-Kahan-Lanczos factorization ($mode)" for mode in + (:vector, :inplace, + :outplace, :mixed) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) v = A * rand(T, (n,)) # ensure v is in column space of A - iter = GKLIterator(A, v, orth) + iter = GKLIterator(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), orth) verbosity = 1 fact = @constinferred initialize(iter; verbosity=verbosity) while length(fact) < n @@ -25,8 +17,8 @@ end verbosity = 0 end - U = stack(basis(fact, :U)) - V = stack(basis(fact, :V)) + U = stack(unwrapvec, basis(fact, :U)) + V = stack(unwrapvec, basis(fact, :V)) B = rayleighquotient(fact) @test normres(fact) < 10 * n * eps(real(T)) @test U' * U ≈ I @@ -39,65 +31,17 @@ end @test rayleighquotient(last(states)) ≈ B end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) - v = MinimalVec{IP}(A * rand(T, (n,))) # ensure v is in column space of A - iter = GKLIterator(wrapop(A), v) - verbosity = 1 - fact = @constinferred initialize(iter; verbosity=verbosity) - while length(fact) < n - @constinferred expand!(iter, fact; verbosity=verbosity) - verbosity = 0 - end - - U = stack(unwrap, basis(fact, :U)) - V = stack(unwrap, basis(fact, :V)) - B = rayleighquotient(fact) - @test normres(fact) < 10 * n * eps(real(T)) - @test V' * V ≈ I - @test U' * U ≈ I - @test A * V ≈ U * B - @test A' * U ≈ V * B' - - @constinferred initialize!(iter, deepcopy(fact); verbosity=1) - states = collect(Iterators.take(iter, n)) # collect tests size and eltype? - @test rayleighquotient(last(states)) ≈ B - end - - @testset "MixedVec" begin - T = ComplexF64 - A = rand(T, (n, n)) - v = MinimalVec{false}(A * rand(T, (n,))) # ensure v is in column space of A - - iter = GKLIterator(wrapop2(A), v) - verbosity = 1 - fact = @constinferred initialize(iter; verbosity=verbosity) - while length(fact) < n - @constinferred expand!(iter, fact; verbosity=verbosity) - verbosity = 0 - end - - U = stack(unwrap, basis(fact, :U)) - V = stack(unwrap, basis(fact, :V)) - B = rayleighquotient(fact) - @test normres(fact) < 10 * n * eps(real(T)) - @test V' * V ≈ I - @test U' * U ≈ I - @test A * V ≈ U * B - @test A' * U ≈ V * B' - - @constinferred initialize!(iter, deepcopy(fact); verbosity=1) - states = collect(Iterators.take(iter, n)) # collect tests size and eltype? - @test rayleighquotient(last(states)) ≈ B - end end # Test incomplete Golub-Kahan-Lanczos factorization -@testset "Incomplete Golub-Kahan-Lanczos factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int}) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Incomplete Golub-Kahan-Lanczos factorization ($mode)" for mode in + (:vector, :inplace, + :outplace, :mixed) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths if T == Complex{Int} A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N)) v = rand(-100:100, (N,)) @@ -105,15 +49,16 @@ end A = rand(T, (N, N)) v = rand(T, (N,)) end - iter = @constinferred GKLIterator(A, v, orth) + iter = @constinferred GKLIterator(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), + orth) krylovdim = 3 * n fact = @constinferred initialize(iter) while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim @constinferred expand!(iter, fact) - - Ũ, Ṽ, B, r, β, e = fact - U = stack(Ũ) - V = stack(Ṽ) + Ũ, Ṽ, B, r̃, β, e = fact + U = stack(unwrapvec, Ũ) + V = stack(unwrapvec, Ṽ) + r = unwrapvec(r̃) @test U' * U ≈ I @test V' * V ≈ I @test norm(r) ≈ β @@ -122,10 +67,10 @@ end end fact = @constinferred shrink!(fact, div(n, 2)) - U = stack(@constinferred basis(fact, :U)) - V = stack(@constinferred basis(fact, :V)) + U = stack(unwrapvec, @constinferred basis(fact, :U)) + V = stack(unwrapvec, @constinferred basis(fact, :V)) B = @constinferred rayleighquotient(fact) - r = @constinferred residual(fact) + r = unwrapvec(@constinferred residual(fact)) β = @constinferred normres(fact) e = @constinferred rayleighextension(fact) @test U' * U ≈ I @@ -135,72 +80,4 @@ end @test A' * U ≈ V * B' end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) - v = MinimalVec{IP}(rand(T, (N,))) - iter = @constinferred GKLIterator(wrapop(A), v) - krylovdim = 3 * n - fact = @constinferred initialize(iter) - while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim - @constinferred expand!(iter, fact) - Ũ, Ṽ, B, r̃, β, e = fact - U = stack(unwrap, Ũ) - V = stack(unwrap, Ṽ) - r = unwrap(r̃) - @test U' * U ≈ I - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ U * B + r * e' - @test A' * U ≈ V * B' - end - - fact = @constinferred shrink!(fact, div(n, 2)) - U = stack(unwrap, @constinferred basis(fact, :U)) - V = stack(unwrap, @constinferred basis(fact, :V)) - B = @constinferred rayleighquotient(fact) - r = unwrap(@constinferred residual(fact)) - β = @constinferred normres(fact) - e = @constinferred rayleighextension(fact) - @test U' * U ≈ I - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ U * B + r * e' - @test A' * U ≈ V * B' - end - - @testset "MixedVec" begin - T = ComplexF64 - A = rand(T, (N, N)) - v = MinimalVec{false}(rand(T, (N,))) - iter = @constinferred GKLIterator(wrapop2(A), v) - krylovdim = 3 * n - fact = @constinferred initialize(iter) - while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim - @constinferred expand!(iter, fact) - Ũ, Ṽ, B, r̃, β, e = fact - U = stack(unwrap, Ũ) - V = stack(unwrap, Ṽ) - r = unwrap(r̃) - @test U' * U ≈ I - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ U * B + r * e' - @test A' * U ≈ V * B' - end - - fact = @constinferred shrink!(fact, div(n, 2)) - U = stack(unwrap, @constinferred basis(fact, :U)) - V = stack(unwrap, @constinferred basis(fact, :V)) - B = @constinferred rayleighquotient(fact) - r = unwrap(@constinferred residual(fact)) - β = @constinferred normres(fact) - e = @constinferred rayleighextension(fact) - @test U' * U ≈ I - @test V' * V ≈ I - @test norm(r) ≈ β - @test A * V ≈ U * B + r * e' - @test A' * U ≈ V * B' - end end diff --git a/test/linsolve.jl b/test/linsolve.jl index 2aa78b5..2246011 100644 --- a/test/linsolve.jl +++ b/test/linsolve.jl @@ -1,17 +1,19 @@ # Test CG complete -@testset "CG small problem" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) +@testset "CG small problem ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + @testset for T in scalartypes A = rand(T, (n, n)) A = sqrt(A * A') b = rand(T, n) alg = CG(; maxiter=2n, tol=tolerance(T) * norm(b), verbosity=2) # because of loss of orthogonality, we choose maxiter = 2n - x, info = @constinferred linsolve(A, b; + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)); ishermitian=true, isposdef=true, maxiter=2n, krylovdim=1, rtol=tolerance(T), verbosity=1) @test info.converged > 0 - @test b ≈ A * x - x, info = @constinferred linsolve(A, b, x; + @test unwrapvec(b) ≈ A * unwrapvec(x) + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), x; ishermitian=true, isposdef=true, maxiter=2n, krylovdim=1, rtol=tolerance(T)) @test info.numops == 1 @@ -20,84 +22,51 @@ A = sqrt(A * A') α₀ = rand(real(T)) + 1 α₁ = rand(real(T)) - x, info = @constinferred linsolve(A, b, zerovector(b), alg, α₀, α₁) - @test b ≈ (α₀ * I + α₁ * A) * x - @test info.converged > 0 - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) - A += A' - b = MinimalVec{IP}(rand(T, (n,))) - alg = CG(; maxiter=2n, tol=tolerance(T) * norm(b), verbosity=2) # because of loss of orthogonality, we choose maxiter = 2n - x, info = @constinferred linsolve(wrapop(A), b; ishermitian=true, isposdef=true, - maxiter=2n, krylovdim=1, rtol=tolerance(T), - verbosity=1) - @test info.converged > 0 - @test unwrap(b) ≈ A * unwrap(x) - x, info = @constinferred linsolve(wrapop(A), b, x; ishermitian=true, isposdef=true, - maxiter=2n, krylovdim=1, rtol=tolerance(T)) - @test info.numops == 1 - - A = rand(T, (n, n)) - A += A' - α₀ = rand(real(T)) + 1 - α₁ = rand(real(T)) - x, info = @constinferred linsolve(wrapop(A), b, zerovector(b), alg, α₀, α₁) - @test unwrap(b) ≈ (α₀ * I + α₁ * A) * unwrap(x) + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), + wrapvec(zerovector(b), Val(mode)), alg, α₀, α₁) + @test unwrapvec(b) ≈ (α₀ * I + α₁ * A) * unwrapvec(x) @test info.converged > 0 end end # Test CG complete -@testset "CG large problem" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) +@testset "CG large problem ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + @testset for T in scalartypes A = rand(T, (N, N)) A = sqrt(sqrt(A * A')) / N b = rand(T, N) - x, info = @constinferred linsolve(A, b; - isposdef=true, maxiter=1, krylovdim=N, - rtol=tolerance(T)) - @test b ≈ A * x + info.residual - - α₀ = rand(real(T)) + 1 - α₁ = rand(real(T)) - x, info = @constinferred linsolve(A, b, α₀, α₁; + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)); isposdef=true, maxiter=1, krylovdim=N, rtol=tolerance(T)) - @test b ≈ (α₀ * I + α₁ * A) * x + info.residual - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) - A += A' - b = MinimalVec{IP}(rand(T, (N,))) - x, info = @constinferred linsolve(wrapop(A), b; isposdef=true, maxiter=1, - krylovdim=N, rtol=tolerance(T)) - @test unwrap(b) ≈ A * unwrap(x) + unwrap(info.residual) + @test unwrapvec(b) ≈ A * unwrapvec(x) + unwrapvec(info.residual) α₀ = rand(real(T)) + 1 α₁ = rand(real(T)) - x, info = @constinferred linsolve(wrapop(A), b, zerovector(b), α₀, α₁; + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), + α₀, α₁; isposdef=true, maxiter=1, krylovdim=N, rtol=tolerance(T)) - @test unwrap(b) ≈ (α₀ * I + α₁ * A) * unwrap(x) + unwrap(info.residual) + @test unwrapvec(b) ≈ (α₀ * I + α₁ * A) * unwrapvec(x) + unwrapvec(info.residual) end end # Test GMRES complete -@testset "GMRES full factorization" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) +@testset "GMRES full factorization ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + @testset for T in scalartypes A = rand(T, (n, n)) .- one(T) / 2 b = rand(T, n) alg = GMRES(; krylovdim=n, maxiter=2, tol=tolerance(T) * norm(b), verbosity=2) - x, info = @constinferred linsolve(A, b; krylovdim=n, maxiter=2, + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)); + krylovdim=n, maxiter=2, rtol=tolerance(T), verbosity=1) @test info.converged > 0 - @test b ≈ A * x - x, info = @constinferred linsolve(A, b, x; krylovdim=n, maxiter=2, + @test unwrapvec(b) ≈ A * unwrapvec(x) + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), x; + krylovdim=n, maxiter=2, rtol=tolerance(T)) @test info.numops == 1 @@ -105,81 +74,50 @@ end α₀ = rand(T) α₁ = -rand(T) x, info = @constinferred(linsolve(A, b, zerovector(b), alg, α₀, α₁)) - @test b ≈ (α₀ * I + α₁ * A) * x - @test info.converged > 0 - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) - b = MinimalVec{IP}(rand(T, (n,))) - alg = GMRES(; krylovdim=n, maxiter=2, tol=tolerance(T) * norm(b), verbosity=2) - x, info = @constinferred linsolve(wrapop(A), b; krylovdim=n, maxiter=2, - rtol=tolerance(T), verbosity=1) - @test info.converged > 0 - @test unwrap(b) ≈ A * unwrap(x) - x, info = @constinferred linsolve(wrapop(A), b, x; krylovdim=n, maxiter=2, - rtol=tolerance(T)) - @test info.numops == 1 - - A = rand(T, (n, n)) - α₀ = rand(T) - α₁ = -rand(T) - x, info = @constinferred(linsolve(wrapop(A), b, zerovector(b), alg, α₀, α₁)) - @test unwrap(b) ≈ (α₀ * I + α₁ * A) * unwrap(x) + @test unwrapvec(b) ≈ (α₀ * I + α₁ * A) * unwrapvec(x) @test info.converged > 0 end end # Test GMRES with restart -@testset "GMRES with restarts" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) +@testset "GMRES with restarts ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + @testset for T in scalartypes A = rand(T, (N, N)) .- one(T) / 2 A = I - T(9 / 10) * A / maximum(abs, eigvals(A)) b = rand(T, N) - x, info = @constinferred linsolve(A, b; krylovdim=3 * n, - maxiter=50, rtol=tolerance(T)) - @test b ≈ A * x + info.residual - - A = rand(T, (N, N)) .- one(T) / 2 - α₀ = maximum(abs, eigvals(A)) - α₁ = -rand(T) - α₁ *= T(9) / T(10) / abs(α₁) - x, info = @constinferred linsolve(A, b, α₀, α₁; krylovdim=3 * n, - maxiter=50, rtol=tolerance(T)) - @test b ≈ (α₀ * I + α₁ * A) * x + info.residual - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - b = MinimalVec{IP}(rand(T, (N,))) - x, info = @constinferred linsolve(wrapop(A), b; krylovdim=3 * n, + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)); + krylovdim=3 * n, maxiter=50, rtol=tolerance(T)) - @test unwrap(b) ≈ A * unwrap(x) + unwrap(info.residual) + @test unwrapvec(b) ≈ A * unwrapvec(x) + unwrapvec(info.residual) A = rand(T, (N, N)) .- one(T) / 2 α₀ = maximum(abs, eigvals(A)) α₁ = -rand(T) α₁ *= T(9) / T(10) / abs(α₁) - x, info = @constinferred linsolve(wrapop(A), b, zerovector(b), α₀, α₁; - krylovdim=3 * n, + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), α₀, + α₁; krylovdim=3 * n, maxiter=50, rtol=tolerance(T)) - @test unwrap(b) ≈ (α₀ * I + α₁ * A) * unwrap(x) + unwrap(info.residual) + @test unwrapvec(b) ≈ (α₀ * I + α₁ * A) * unwrapvec(x) + unwrapvec(info.residual) end end # Test BICGStab -@testset "BiCGStab" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) +@testset "BiCGStab ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + @testset for T in scalartypes A = rand(T, (n, n)) .- one(T) / 2 A = I - T(9 / 10) * A / maximum(abs, eigvals(A)) b = rand(T, n) alg = BiCGStab(; maxiter=4n, tol=tolerance(T) * norm(b), verbosity=2) - x, info = @constinferred linsolve(A, b, zerovector(b), alg) + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), + wrapvec(zerovector(b), Val(mode)), alg) @test info.converged > 0 - @test b ≈ A * x - x, info = @constinferred linsolve(A, b, x, alg) + @test unwrapvec(b) ≈ A * unwrapvec(x) + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), x, + alg) @test info.numops == 1 A = rand(T, (N, N)) .- one(T) / 2 @@ -188,37 +126,14 @@ end α₁ = -rand(T) α₁ *= T(9) / T(10) / abs(α₁) alg = BiCGStab(; maxiter=2, tol=tolerance(T) * norm(b), verbosity=1) - x, info = @constinferred linsolve(A, b, zerovector(b), alg, α₀, + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), + wrapvec(zerovector(b), Val(mode)), alg, α₀, α₁) - @test b ≈ (α₀ * I + α₁ * A) * x + info.residual - alg = BiCGStab(; maxiter=10 * N, tol=tolerance(T) * norm(b), verbosity=0) - x, info = @constinferred linsolve(A, b, x, alg, α₀, α₁) - @test info.converged > 0 - @test b ≈ (α₀ * I + α₁ * A) * x - end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - b = MinimalVec{IP}(rand(T, (n,))) - alg = BiCGStab(; maxiter=4n, tol=tolerance(T) * norm(b), verbosity=2) - x, info = @constinferred linsolve(wrapop(A), b, zerovector(b), alg) - @test info.converged > 0 - @test unwrap(b) ≈ A * unwrap(x) - x, info = @constinferred linsolve(wrapop(A), b, x, alg) - @test info.numops == 1 - - A = rand(T, (N, N)) .- one(T) / 2 - b = MinimalVec{IP}(rand(T, (N,))) - α₀ = maximum(abs, eigvals(A)) - α₁ = -rand(T) - α₁ *= T(9) / T(10) / abs(α₁) - alg = BiCGStab(; maxiter=2, tol=tolerance(T) * norm(b), verbosity=1) - x, info = @constinferred linsolve(wrapop(A), b, zerovector(b), alg, α₀, α₁) - @test unwrap(b) ≈ (α₀ * I + α₁ * A) * unwrap(x) + unwrap(info.residual) + @test unwrapvec(b) ≈ (α₀ * I + α₁ * A) * unwrapvec(x) + unwrapvec(info.residual) alg = BiCGStab(; maxiter=10 * N, tol=tolerance(T) * norm(b), verbosity=0) - x, info = @constinferred linsolve(wrapop(A), b, x, alg, α₀, α₁) + x, info = @constinferred linsolve(wrapop(A, Val(mode)), wrapvec(b, Val(mode)), x, + alg, α₀, α₁) @test info.converged > 0 - @test unwrap(b) ≈ (α₀ * I + α₁ * A) * unwrap(x) + @test unwrapvec(b) ≈ (α₀ * I + α₁ * A) * unwrapvec(x) end end diff --git a/test/schursolve.jl b/test/schursolve.jl index 026794f..aebd418 100644 --- a/test/schursolve.jl +++ b/test/schursolve.jl @@ -1,20 +1,26 @@ -@testset "Arnoldi - schursolve full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Arnoldi - schursolve full ($mode)" for mode in (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) .- one(T) / 2 v = rand(T, (n,)) alg = Arnoldi(; orth=orth, krylovdim=n, maxiter=1, tol=tolerance(T)) n1 = div(n, 2) - T1, V1, D1, info1 = @constinferred schursolve(A, v, n1, :SR, alg) + T1, V1, D1, info1 = @constinferred schursolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), n1, :SR, + alg) n2 = n - n1 - T2, V2, D2, info2 = schursolve(A, v, n2, :LR, alg) + T2, V2, D2, info2 = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n2, + :LR, alg) D = sort(sort(eigvals(A); by=imag, rev=true); alg=MergeSort, by=real) D2′ = sort(sort(D2; by=imag, rev=true); alg=MergeSort, by=real) @test vcat(D1[1:n1], D2′[(end - n2 + 1):end]) ≈ D - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1' * U1 ≈ I @test U2' * U2 ≈ I @@ -22,14 +28,16 @@ @test A * U2 ≈ U2 * T2 if T <: Complex - T1, V1, D1, info = schursolve(A, v, n1, :SI, alg) - T2, V2, D2, info = schursolve(A, v, n2, :LI, alg) + T1, V1, D1, info = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), + n1, :SI, alg) + T2, V2, D2, info = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), + n2, :LI, alg) D = sort(eigvals(A); by=imag) @test vcat(D1[1:n1], reverse(D2[1:n2])) ≈ D - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1' * U1 ≈ I @test U2' * U2 ≈ I @@ -38,54 +46,25 @@ end end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) .- one(T) / 2 - v = MinimalVec{IP}(rand(T, (n,))) - alg = Arnoldi(; orth=cgs2, krylovdim=n, maxiter=1, tol=tolerance(T)) - n1 = div(n, 2) - T1, V1, D1, info1 = @constinferred schursolve(wrapop(A), v, n1, :SR, alg) - n2 = n - n1 - T2, V2, D2, info2 = schursolve(wrapop(A), v, n2, :LR, alg) - D = sort(sort(eigvals(A); by=imag, rev=true); alg=MergeSort, by=real) - D2′ = sort(sort(D2; by=imag, rev=true); alg=MergeSort, by=real) - - @test vcat(D1[1:n1], D2′[(end - n2 + 1):end]) ≈ D - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1' * U1 ≈ I - @test U2' * U2 ≈ I - - @test A * U1 ≈ U1 * T1 - @test A * U2 ≈ U2 * T2 - - T1, V1, D1, info = schursolve(wrapop(A), v, n1, :SI, alg) - T2, V2, D2, info = schursolve(wrapop(A), v, n2, :LI, alg) - D = sort(eigvals(A); by=imag) - - @test vcat(D1[1:n1], reverse(D2[1:n2])) ≈ D - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1' * U1 ≈ I - @test U2' * U2 ≈ I - - @test A * U1 ≈ U1 * T1 - @test A * U2 ≈ U2 * T2 - end end -@testset "Arnoldi - schursolve iteratively" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "Arnoldi - schursolve iteratively ($mode)" for mode in + (:vector, :inplace, :outplace) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (N, N)) .- one(T) / 2 v = rand(T, (N,)) alg = Arnoldi(; orth=orth, krylovdim=3 * n, maxiter=10, tol=tolerance(T)) - T1, V1, D1, info1 = @constinferred schursolve(A, v, n, :SR, alg) - T2, V2, D2, info2 = schursolve(A, v, n, :LR, alg) - T3, V3, D3, info3 = schursolve(A, v, n, :LM, alg) + T1, V1, D1, info1 = @constinferred schursolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), n, :SR, + alg) + T2, V2, D2, info2 = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, + :LR, alg) + T3, V3, D3, info3 = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), n, + :LM, alg) D = sort(eigvals(A); by=imag, rev=true) l1 = info1.converged @@ -95,23 +74,25 @@ end @test D2[1:l2] ≊ sort(D; alg=MergeSort, by=real, rev=true)[1:l2] @test D3[1:l3] ≊ sort(D; alg=MergeSort, by=abs, rev=true)[1:l3] - U1 = stack(V1) - U2 = stack(V2) - U3 = stack(V3) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) + U3 = stack(unwrapvec, V3) @test U1' * U1 ≈ one(U1' * U1) @test U2' * U2 ≈ one(U2' * U2) @test U3' * U3 ≈ one(U3' * U3) - R1 = stack(info1.residual) - R2 = stack(info2.residual) - R3 = stack(info3.residual) + R1 = stack(unwrapvec, info1.residual) + R2 = stack(unwrapvec, info2.residual) + R3 = stack(unwrapvec, info3.residual) @test A * U1 ≈ U1 * T1 + R1 @test A * U2 ≈ U2 * T2 + R2 @test A * U3 ≈ U3 * T3 + R3 if T <: Complex - T1, V1, D1, info1 = schursolve(A, v, n, :SI, alg) - T2, V2, D2, info2 = schursolve(A, v, n, :LI, alg) + T1, V1, D1, info1 = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), + n, :SI, alg) + T2, V2, D2, info2 = schursolve(wrapop(A, Val(mode)), wrapvec(v, Val(mode)), + n, :LI, alg) D = eigvals(A) l1 = info1.converged @@ -119,67 +100,16 @@ end @test D1[1:l1] ≊ sort(D; by=imag)[1:l1] @test D2[1:l2] ≊ sort(D; by=imag, rev=true)[1:l2] - U1 = stack(V1) - U2 = stack(V2) + U1 = stack(unwrapvec, V1) + U2 = stack(unwrapvec, V2) @test U1[:, 1:l1]' * U1[:, 1:l1] ≈ I @test U2[:, 1:l2]' * U2[:, 1:l2] ≈ I - R1 = stack(info1.residual) - R2 = stack(info2.residual) + R1 = stack(unwrapvec, info1.residual) + R2 = stack(unwrapvec, info2.residual) @test A * U1 ≈ U1 * T1 + R1 @test A * U2 ≈ U2 * T2 + R2 end end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (N, N)) .- one(T) / 2 - v = MinimalVec{IP}(rand(T, (N,))) - alg = Arnoldi(; orth=cgs2, krylovdim=3 * n, maxiter=10, tol=tolerance(T)) - T1, V1, D1, info1 = @constinferred schursolve(wrapop(A), v, n, :SR, alg) - T2, V2, D2, info2 = schursolve(wrapop(A), v, n, :LR, alg) - T3, V3, D3, info3 = schursolve(wrapop(A), v, n, :LM, alg) - D = sort(eigvals(A); by=imag, rev=true) - - l1 = info1.converged - l2 = info2.converged - l3 = info3.converged - @test D1[1:l1] ≊ sort(D; alg=MergeSort, by=real)[1:l1] - @test D2[1:l2] ≊ sort(D; alg=MergeSort, by=real, rev=true)[1:l2] - @test D3[1:l3] ≊ sort(D; alg=MergeSort, by=abs, rev=true)[1:l3] - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - U3 = stack(unwrap, V3) - @test U1' * U1 ≈ one(U1' * U1) - @test U2' * U2 ≈ one(U2' * U2) - @test U3' * U3 ≈ one(U3' * U3) - - R1 = stack(unwrap, info1.residual) - R2 = stack(unwrap, info2.residual) - R3 = stack(unwrap, info3.residual) - @test A * U1 ≈ U1 * T1 + R1 - @test A * U2 ≈ U2 * T2 + R2 - @test A * U3 ≈ U3 * T3 + R3 - - T1, V1, D1, info1 = schursolve(wrapop(A), v, n, :SI, alg) - T2, V2, D2, info2 = schursolve(wrapop(A), v, n, :LI, alg) - D = eigvals(A) - - l1 = info1.converged - l2 = info2.converged - @test D1[1:l1] ≊ sort(D; by=imag)[1:l1] - @test D2[1:l2] ≊ sort(D; by=imag, rev=true)[1:l2] - - U1 = stack(unwrap, V1) - U2 = stack(unwrap, V2) - @test U1[:, 1:l1]' * U1[:, 1:l1] ≈ I - @test U2[:, 1:l2]' * U2[:, 1:l2] ≈ I - - R1 = stack(unwrap, info1.residual) - R2 = stack(unwrap, info2.residual) - @test A * U1 ≈ U1 * T1 + R1 - @test A * U2 ≈ U2 * T2 + R2 - end end diff --git a/test/svdsolve.jl b/test/svdsolve.jl index f831b55..bad5157 100644 --- a/test/svdsolve.jl +++ b/test/svdsolve.jl @@ -1,118 +1,52 @@ -@testset "GKL - svdsolve full" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "GKL - svdsolve full ($mode)" for mode in (:vector, :inplace, :outplace, :mixed) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (n, n)) alg = GKL(; orth=orth, krylovdim=2 * n, maxiter=1, tol=tolerance(T)) - S, lvecs, rvecs, info = @constinferred svdsolve(A, A[:, 1], n, :LR, alg) + S, lvecs, rvecs, info = @constinferred svdsolve(wrapop(A, Val(mode)), + wrapvec(A[:, 1], Val(mode)), n, + :LR, alg) @test S ≈ svdvals(A) - U = stack(lvecs) - V = stack(rvecs) + U = stack(unwrapvec, lvecs) + V = stack(unwrapvec, rvecs) @test U' * U ≈ I @test V' * V ≈ I @test A * V ≈ U * Diagonal(S) end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (n, n)) - v = MinimalVec{IP}(rand(T, (n,))) - alg = GKL(; krylovdim=2 * n, maxiter=1, tol=tolerance(T)) - S, lvecs, rvecs, info = @constinferred svdsolve(wrapop(A), v, n, :LR, alg) - - @test S ≈ svdvals(A) - - U = stack(unwrap, lvecs) - V = stack(unwrap, rvecs) - @test U' * U ≈ I - @test V' * V ≈ I - @test A * V ≈ U * Diagonal(S) - end - - @testset "MixedVec" begin - T = ComplexF64 - A = rand(T, (n, n)) - v = MinimalVec{false}(rand(T, (n,))) - alg = GKL(; krylovdim=2 * n, maxiter=1, tol=tolerance(T)) - S, lvecs, rvecs, info = @constinferred svdsolve(wrapop2(A), v, n, :LR, alg) - - @test S ≈ svdvals(A) - - U = stack(unwrap, lvecs) - V = stack(unwrap, rvecs) - @test U' * U ≈ I - @test V' * V ≈ I - @test A * V ≈ U * Diagonal(S) - end end -@testset "GKL - svdsolve iteratively" begin - @testset for T in (Float32, Float64, ComplexF32, ComplexF64) - @testset for orth in (cgs2, mgs2, cgsr, mgsr) +@testset "GKL - svdsolve iteratively ($mode)" for mode in + (:vector, :inplace, :outplace, :mixed) + scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : + (ComplexF64,) + orths = mode === :vector ? (cgs2, mgs2, cgsr, mgsr) : (mgsr,) + @testset for T in scalartypes + @testset for orth in orths A = rand(T, (2 * N, N)) v = rand(T, (2 * N,)) n₁ = div(n, 2) alg = GKL(; orth=orth, krylovdim=n, maxiter=10, tol=tolerance(T), eager=true) - S, lvecs, rvecs, info = @constinferred svdsolve(A, v, n₁, :LR, - alg) + S, lvecs, rvecs, info = @constinferred svdsolve(wrapop(A, Val(mode)), + wrapvec(v, Val(mode)), + n₁, :LR, alg) l = info.converged @test S[1:l] ≈ svdvals(A)[1:l] - U = stack(lvecs) - V = stack(rvecs) + U = stack(unwrapvec, lvecs) + V = stack(unwrapvec, rvecs) @test U[:, 1:l]' * U[:, 1:l] ≈ I @test V[:, 1:l]' * V[:, 1:l] ≈ I - R = stack(info.residual) + R = stack(unwrapvec, info.residual) @test A' * U ≈ V * Diagonal(S) @test A * V ≈ U * Diagonal(S) + R end end - - @testset "MinimalVec{$IP}" for IP in (true, false) - T = ComplexF64 - A = rand(T, (2 * N, N)) - v = MinimalVec{IP}(rand(T, (2 * N,))) - n₁ = div(n, 2) - alg = GKL(; krylovdim=n, maxiter=10, tol=tolerance(T), eager=true) - S, lvecs, rvecs, info = @constinferred svdsolve(wrapop(A), v, n₁, :LR, - alg) - - l = info.converged - @test S[1:l] ≈ svdvals(A)[1:l] - - U = stack(unwrap, lvecs) - V = stack(unwrap, rvecs) - @test U[:, 1:l]' * U[:, 1:l] ≈ I - @test V[:, 1:l]' * V[:, 1:l] ≈ I - - R = stack(unwrap, info.residual) - @test A' * U ≈ V * Diagonal(S) - @test A * V ≈ U * Diagonal(S) + R - end - - @testset "MixedVec" begin - T = ComplexF64 - A = rand(T, (2 * N, N)) - v = MinimalVec{false}(rand(T, (2 * N,))) - n₁ = div(n, 2) - alg = GKL(; krylovdim=n, maxiter=10, tol=tolerance(T), eager=true) - S, lvecs, rvecs, info = @constinferred svdsolve(wrapop2(A), v, n₁, :LR, - alg) - - l = info.converged - @test S[1:l] ≈ svdvals(A)[1:l] - - U = stack(unwrap, lvecs) - V = stack(unwrap, rvecs) - @test U[:, 1:l]' * U[:, 1:l] ≈ I - @test V[:, 1:l]' * V[:, 1:l] ≈ I - - R = stack(unwrap, info.residual) - @test A' * U ≈ V * Diagonal(S) - @test A * V ≈ U * Diagonal(S) + R - end end diff --git a/test/testsetup.jl b/test/testsetup.jl index ad59bb2..8d4bb91 100644 --- a/test/testsetup.jl +++ b/test/testsetup.jl @@ -1,6 +1,7 @@ module TestSetup -export tolerance, ≊, MinimalVec, unwrap, isinplace, stack +export tolerance, ≊, MinimalVec, isinplace, stack +export wrapop, wrapvec, unwrapvec import VectorInterface as VI using VectorInterface @@ -45,8 +46,6 @@ const OutplaceVec{V} = MinimalVec{false,V} isinplace(::Type{MinimalVec{IP,V}}) where {V,IP} = IP isinplace(v::MinimalVec) = isinplace(typeof(v)) -unwrap(v::MinimalVec) = v.vec - VI.scalartype(::Type{<:MinimalVec{IP,V}}) where {IP,V} = scalartype(V) function VI.zerovector(v::MinimalVec, S::Type{<:Number}) @@ -91,8 +90,44 @@ end VI.inner(x::MinimalVec, y::MinimalVec) = inner(x.vec, y.vec) VI.norm(x::MinimalVec) = LinearAlgebra.norm(x.vec) +# Wrappers +# -------- +# dispatch on val is necessary for type stability + +function wrapvec(v, ::Val{mode}) where {mode} + return mode === :vector ? v : + mode === :inplace ? MinimalVec{true}(v) : + mode === :outplace ? MinimalVec{false}(v) : + mode === :mixed ? MinimalVec{false}(v) : + throw(ArgumentError("invalid mode ($mode)")) +end +function wrapvec2(v, ::Val{mode}) where {mode} + return mode === :mixed ? MinimalVec{true}(v) : wrapvec(v, mode) +end + +unwrapvec(v::MinimalVec) = v.vec +unwrapvec(v) = v + +function wrapop(A, ::Val{mode}) where {mode} + if mode === :vector + return A + elseif mode === :inplace || mode === :outplace + return function (v, flag=Val(false)) + if flag === Val(true) + return wrapvec(A' * unwrapvec(v), Val(mode)) + else + return wrapvec(A * unwrapvec(v), Val(mode)) + end + end + elseif mode === :mixed + return (x -> wrapvec(A * unwrapvec(x), Val(mode)), + y -> wrapvec2(A' * unwrapvec(y), Val(mode))) + else + throw(ArgumentError("invalid mode ($mode)")) + end +end + if VERSION < v"1.9" - stack(itr) = reduce(hcat, itr) stack(f, itr) = mapreduce(f, hcat, itr) end