From e86af5b58247007c2191dce57c3a4f66fb2224d1 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sun, 8 Oct 2017 16:56:41 +0200 Subject: [PATCH] reset GLOBAL_RNG state in/out at-testset for reproducibility This is a follow-up on the recently introduced `guardsrand` functionality, first suggested at https://github.com/JuliaLang/julia/pull/16940#issuecomment-226152878. Each `testset` block is now implicitly wrapped in a `guardsrand` block, which is more user-friendly and more systematic (only few users know about `guardsrand`). These are essentially two new features: 1) "in": in base, tests are run with the global RNG randomly seeded, but the seed is printed in case of failure to allow reproducing the failure; but even if the failure occurs after many tests, when they alter the global RNG state, one has to re-run the whole file, which can be time-consuming; with this change, at the beginning of each `testset`, the global RNG is re-seeded with its own seed: this allows to re-run only the failing `testset`, which can be done easily in the REPL (the seeding occurs for each loop in a `testset for`; this also allows to re-arrange `testset`s in arbitrary order w.r.t. the global RNG; 2) "out": a `testset` leaves no tracks of its use of `srand` or `rand` (this "feature" should be less and less needed/useful with the generalization of the use of `testset` blocks). Example: ``` @testset begin srand(123) rand() @testset for T in (Int, Float64) # here is an implicit srand(123), by 1) rand() end rand() # this value will not be affected if the sub-`testset` block # above is removed, or if another rand() call is added therein, by 2) end ``` Note that guardsrand can't be used directly, as then the testset's body is wrapped in a function, which causes problem with overwriting loop variable ("outer"), using `using`, defining new methods, etc. So we need to duplicate guardsrand's logic. --- stdlib/IterativeEigensolvers/test/runtests.jl | 188 +++--- stdlib/Test/src/Test.jl | 14 +- stdlib/Test/test/runtests.jl | 25 +- test/linalg/lapack.jl | 94 ++- test/linalg/lu.jl | 1 + test/linalg/tridiag.jl | 560 +++++++++--------- test/random.jl | 4 +- test/sorting.jl | 118 ++-- test/sparse/sparse.jl | 74 ++- 9 files changed, 551 insertions(+), 527 deletions(-) diff --git a/stdlib/IterativeEigensolvers/test/runtests.jl b/stdlib/IterativeEigensolvers/test/runtests.jl index 4d59d32c513ad..f77d952219ae9 100644 --- a/stdlib/IterativeEigensolvers/test/runtests.jl +++ b/stdlib/IterativeEigensolvers/test/runtests.jl @@ -4,118 +4,118 @@ using IterativeEigensolvers using Test @testset "eigs" begin - guardsrand(1234) do - n = 10 - areal = sprandn(n,n,0.4) - breal = sprandn(n,n,0.4) - acmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4)) - bcmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4)) + srand(1234) + n = 10 + areal = sprandn(n,n,0.4) + breal = sprandn(n,n,0.4) + acmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4)) + bcmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4)) - testtol = 1e-6 + testtol = 1e-6 - @testset for elty in (Float64, ComplexF64) - if elty == ComplexF32 || elty == ComplexF64 - a = acmplx - b = bcmplx - else - a = areal - b = breal - end - a_evs = eigvals(Array(a)) - a = convert(SparseMatrixCSC{elty}, a) - asym = a' + a # symmetric indefinite - apd = a'*a # symmetric positive-definite + @testset for elty in (Float64, ComplexF64) + if elty == ComplexF32 || elty == ComplexF64 + a = acmplx + b = bcmplx + else + a = areal + b = breal + end + a_evs = eigvals(Array(a)) + a = convert(SparseMatrixCSC{elty}, a) + asym = a' + a # symmetric indefinite + apd = a'*a # symmetric positive-definite - b = convert(SparseMatrixCSC{elty}, b) - bsym = b' + b - bpd = b'*b + b = convert(SparseMatrixCSC{elty}, b) + bsym = b' + b + bpd = b'*b - (d,v) = eigs(a, nev=3) - @test a*v[:,2] ≈ d[2]*v[:,2] - @test norm(v) > testtol # eigenvectors cannot be null vectors - (d,v) = eigs(a, I, nev=3) # test eigs(A, B; kwargs...) - @test a*v[:,2] ≈ d[2]*v[:,2] - @test norm(v) > testtol # eigenvectors cannot be null vectors - @test_logs (:warn,"Use symbols instead of strings for specifying which eigenvalues to compute") eigs(a, which="LM") - @test_logs (:warn,"Adjusting ncv from 1 to 4") eigs(a, ncv=1, nev=2) - @test_logs (:warn,"Adjusting nev from $n to $(n-2)") eigs(a, nev=n) - # (d,v) = eigs(a, b, nev=3, tol=1e-8) # not handled yet - # @test a*v[:,2] ≈ d[2]*b*v[:,2] atol=testtol - # @test norm(v) > testtol # eigenvectors cannot be null vectors - if elty <: Base.LinAlg.BlasComplex - sr_ind = indmin(real.(a_evs)) - (d, v) = eigs(a, nev=1, which=:SR) - @test d[1] ≈ a_evs[sr_ind] - si_ind = indmin(imag.(a_evs)) - (d, v) = eigs(a, nev=1, which=:SI) - @test d[1] ≈ a_evs[si_ind] - lr_ind = indmax(real.(a_evs)) - (d, v) = eigs(a, nev=1, which=:LR) - @test d[1] ≈ a_evs[lr_ind] - li_ind = indmax(imag.(a_evs)) - (d, v) = eigs(a, nev=1, which=:LI) - @test d[1] ≈ a_evs[li_ind] - end + (d,v) = eigs(a, nev=3) + @test a*v[:,2] ≈ d[2]*v[:,2] + @test norm(v) > testtol # eigenvectors cannot be null vectors + (d,v) = eigs(a, I, nev=3) # test eigs(A, B; kwargs...) + @test a*v[:,2] ≈ d[2]*v[:,2] + @test norm(v) > testtol # eigenvectors cannot be null vectors + @test_logs (:warn,"Use symbols instead of strings for specifying which eigenvalues to compute") eigs(a, which="LM") + @test_logs (:warn,"Adjusting ncv from 1 to 4") eigs(a, ncv=1, nev=2) + @test_logs (:warn,"Adjusting nev from $n to $(n-2)") eigs(a, nev=n) + # (d,v) = eigs(a, b, nev=3, tol=1e-8) # not handled yet + # @test a*v[:,2] ≈ d[2]*b*v[:,2] atol=testtol + # @test norm(v) > testtol # eigenvectors cannot be null vectors + if elty <: Base.LinAlg.BlasComplex + sr_ind = indmin(real.(a_evs)) + (d, v) = eigs(a, nev=1, which=:SR) + @test d[1] ≈ a_evs[sr_ind] + si_ind = indmin(imag.(a_evs)) + (d, v) = eigs(a, nev=1, which=:SI) + @test d[1] ≈ a_evs[si_ind] + lr_ind = indmax(real.(a_evs)) + (d, v) = eigs(a, nev=1, which=:LR) + @test d[1] ≈ a_evs[lr_ind] + li_ind = indmax(imag.(a_evs)) + (d, v) = eigs(a, nev=1, which=:LI) + @test d[1] ≈ a_evs[li_ind] + end - (d,v) = eigs(asym, nev=3) - @test asym*v[:,1] ≈ d[1]*v[:,1] - @test eigs(asym; nev=1, sigma=d[3])[1][1] ≈ d[3] - @test norm(v) > testtol # eigenvectors cannot be null vectors + (d,v) = eigs(asym, nev=3) + @test asym*v[:,1] ≈ d[1]*v[:,1] + @test eigs(asym; nev=1, sigma=d[3])[1][1] ≈ d[3] + @test norm(v) > testtol # eigenvectors cannot be null vectors - (d,v) = eigs(apd, nev=3) - @test apd*v[:,3] ≈ d[3]*v[:,3] - @test eigs(apd; nev=1, sigma=d[3])[1][1] ≈ d[3] + (d,v) = eigs(apd, nev=3) + @test apd*v[:,3] ≈ d[3]*v[:,3] + @test eigs(apd; nev=1, sigma=d[3])[1][1] ≈ d[3] - (d,v) = eigs(apd, bpd, nev=3, tol=1e-8) - @test apd*v[:,2] ≈ d[2]*bpd*v[:,2] atol=testtol - @test norm(v) > testtol # eigenvectors cannot be null vectors + (d,v) = eigs(apd, bpd, nev=3, tol=1e-8) + @test apd*v[:,2] ≈ d[2]*bpd*v[:,2] atol=testtol + @test norm(v) > testtol # eigenvectors cannot be null vectors - @testset "(shift-and-)invert mode" begin - (d,v) = eigs(apd, nev=3, sigma=0) - @test apd*v[:,3] ≈ d[3]*v[:,3] - @test norm(v) > testtol # eigenvectors cannot be null vectors + @testset "(shift-and-)invert mode" begin + (d,v) = eigs(apd, nev=3, sigma=0) + @test apd*v[:,3] ≈ d[3]*v[:,3] + @test norm(v) > testtol # eigenvectors cannot be null vectors - (d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8) - @test apd*v[:,1] ≈ d[1]*bpd*v[:,1] atol=testtol - @test norm(v) > testtol # eigenvectors cannot be null vectors - end + (d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8) + @test apd*v[:,1] ≈ d[1]*bpd*v[:,1] atol=testtol + @test norm(v) > testtol # eigenvectors cannot be null vectors + end - @testset "ArgumentErrors" begin - @test_throws ArgumentError eigs(rand(elty,2,2)) - @test_throws ArgumentError eigs(a, nev=-1) - @test_throws ArgumentError eigs(a, which=:Z) - @test_throws ArgumentError eigs(a, which=:BE) - @test_throws DimensionMismatch eigs(a, v0=zeros(elty,n+2)) - @test_throws ArgumentError eigs(a, v0=zeros(Int,n)) - if elty == Float64 - @test_throws ArgumentError eigs(a+a.',which=:SI) - @test_throws ArgumentError eigs(a+a.',which=:LI) - @test_throws ArgumentError eigs(a,sigma=rand(ComplexF32)) - end + @testset "ArgumentErrors" begin + @test_throws ArgumentError eigs(rand(elty,2,2)) + @test_throws ArgumentError eigs(a, nev=-1) + @test_throws ArgumentError eigs(a, which=:Z) + @test_throws ArgumentError eigs(a, which=:BE) + @test_throws DimensionMismatch eigs(a, v0=zeros(elty,n+2)) + @test_throws ArgumentError eigs(a, v0=zeros(Int,n)) + if elty == Float64 + @test_throws ArgumentError eigs(a+a.',which=:SI) + @test_throws ArgumentError eigs(a+a.',which=:LI) + @test_throws ArgumentError eigs(a,sigma=rand(ComplexF32)) end end + end - @testset "Symmetric generalized with singular B" begin - n = 10 - k = 3 - A = randn(n,n); A = A'A - B = randn(n,k); B = B*B' - @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k]) - end + @testset "Symmetric generalized with singular B" begin + srand(124) + n = 10 + k = 3 + A = randn(n,n); A = A'A + B = randn(n,k); B = B*B' + @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k]) end end # Problematic example from #6965A let A6965 = [ - 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 - -1.0 2.0 0.0 0.0 0.0 0.0 0.0 1.0 - -1.0 0.0 3.0 0.0 0.0 0.0 0.0 1.0 - -1.0 0.0 0.0 4.0 0.0 0.0 0.0 1.0 - -1.0 0.0 0.0 0.0 5.0 0.0 0.0 1.0 - -1.0 0.0 0.0 0.0 0.0 6.0 0.0 1.0 - -1.0 0.0 0.0 0.0 0.0 0.0 7.0 1.0 - -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 8.0 - ] + 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 + -1.0 2.0 0.0 0.0 0.0 0.0 0.0 1.0 + -1.0 0.0 3.0 0.0 0.0 0.0 0.0 1.0 + -1.0 0.0 0.0 4.0 0.0 0.0 0.0 1.0 + -1.0 0.0 0.0 0.0 5.0 0.0 0.0 1.0 + -1.0 0.0 0.0 0.0 0.0 6.0 0.0 1.0 + -1.0 0.0 0.0 0.0 0.0 0.0 7.0 1.0 + -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 8.0 +] d, = eigs(A6965,which=:SM,nev=2,ncv=4,tol=eps()) @test d[1] ≈ 2.5346936860350002 @test real(d[2]) ≈ 2.6159972444834976 diff --git a/stdlib/Test/src/Test.jl b/stdlib/Test/src/Test.jl index a848bd85eb52d..35562a8ebcbb0 100644 --- a/stdlib/Test/src/Test.jl +++ b/stdlib/Test/src/Test.jl @@ -964,12 +964,20 @@ function testset_beginend(args, tests, source) # which is needed for backtrace scrubbing to work correctly. while false; end push_testset(ts) + # we reproduce the logic of guardsrand, but this function + # cannot be used as it changes slightly the semantic of @testset, + # by wrapping the body in a function + oldrng = copy(Base.GLOBAL_RNG) try + # GLOBAL_RNG is re-seeded with its own seed to ease reproduce a failed test + srand(Base.GLOBAL_RNG.seed) $(esc(tests)) catch err # something in the test block threw an error. Count that as an # error in this test set record(ts, Error(:nontest_error, :(), err, catch_backtrace(), $(QuoteNode(source)))) + finally + copy!(Base.GLOBAL_RNG, oldrng) end pop_testset() finish(ts) @@ -1030,12 +1038,16 @@ function testset_forloop(args, testloop, source) ts = $(testsettype)($desc; $options...) push_testset(ts) first_iteration = false + oldrng = copy(Base.GLOBAL_RNG) try + srand(Base.GLOBAL_RNG.seed) $(esc(tests)) catch err # Something in the test block threw an error. Count that as an # error in this test set record(ts, Error(:nontest_error, :(), err, catch_backtrace(), $(QuoteNode(source)))) + finally + copy!(Base.GLOBAL_RNG, oldrng) end end quote @@ -1477,7 +1489,7 @@ end "`guardsrand(f, seed)` is equivalent to running `srand(seed); f()` and then restoring the state of the global RNG as it was before." -guardsrand(f::Function, seed::Integer) = guardsrand() do +guardsrand(f::Function, seed::Union{Vector{UInt32},Integer}) = guardsrand() do srand(seed) f() end diff --git a/stdlib/Test/test/runtests.jl b/stdlib/Test/test/runtests.jl index 8b745114e9e21..4a5d8fb76e7c9 100644 --- a/stdlib/Test/test/runtests.jl +++ b/stdlib/Test/test/runtests.jl @@ -287,10 +287,8 @@ end end end @testset "some loops fail" begin - guardsrand(123) do - @testset for i in 1:5 - @test i <= rand(1:10) - end + @testset for i in 1:5 + @test i <= 4 end # should add 3 errors and 3 passing tests @testset for i in 1:6 @@ -739,3 +737,22 @@ end be tested in --depwarn=error mode""" end end + +@testset "@testset preserves GLOBAL_RNG's state, and re-seeds it" begin + # i.e. it behaves as if it was wrapped in a `guardsrand(GLOBAL_RNG.seed)` block + seed = rand(UInt128) + srand(seed) + a = rand() + @testset begin + # global RNG must re-seeded at the beginning of @testset + @test a == rand() + end + @testset for i=1:3 + @test a == rand() + end + # the @testset's above must have no consequence for rand() below + b = rand() + srand(seed) + @test a == rand() + @test b == rand() +end diff --git a/test/linalg/lapack.jl b/test/linalg/lapack.jl index 889841935671a..b23954d09a1ae 100644 --- a/test/linalg/lapack.jl +++ b/test/linalg/lapack.jl @@ -10,25 +10,24 @@ import Base.LinAlg.BlasInt @test_throws ArgumentError Base.LinAlg.LAPACK.chktrans('Z') @testset "syevr" begin - guardsrand(123) do - Ainit = randn(5,5) - @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) - if elty == ComplexF32 || elty == ComplexF64 - A = complex.(Ainit, Ainit) - else - A = Ainit - end - A = convert(Array{elty, 2}, A) - Asym = A'A - vals, Z = LAPACK.syevr!('V', copy(Asym)) - @test Z*(Diagonal(vals)*Z') ≈ Asym - @test all(vals .> 0.0) - @test LAPACK.syevr!('N','V','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[vals .< 1.0] - @test LAPACK.syevr!('N','I','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[4:5] - @test vals ≈ LAPACK.syev!('N','U',copy(Asym)) - - @test_throws DimensionMismatch LAPACK.sygvd!(1,'V','U',copy(Asym),ones(elty,6,6)) + srand(123) + Ainit = randn(5,5) + @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) + if elty == ComplexF32 || elty == ComplexF64 + A = complex.(Ainit, Ainit) + else + A = Ainit end + A = convert(Array{elty, 2}, A) + Asym = A'A + vals, Z = LAPACK.syevr!('V', copy(Asym)) + @test Z*(Diagonal(vals)*Z') ≈ Asym + @test all(vals .> 0.0) + @test LAPACK.syevr!('N','V','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[vals .< 1.0] + @test LAPACK.syevr!('N','I','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[4:5] + @test vals ≈ LAPACK.syev!('N','U',copy(Asym)) + + @test_throws DimensionMismatch LAPACK.sygvd!(1,'V','U',copy(Asym),ones(elty,6,6)) end end @@ -207,12 +206,11 @@ end @testset "gels" begin @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) - guardsrand(913) do - A = rand(elty,10,10) - X = rand(elty,10) - B,Y,z = LAPACK.gels!('N',copy(A),copy(X)) - @test A\X ≈ Y - end + srand(913) + A = rand(elty,10,10) + X = rand(elty,10) + B,Y,z = LAPACK.gels!('N',copy(A),copy(X)) + @test A\X ≈ Y end end @@ -435,36 +433,34 @@ end @testset "sysv" begin @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) - guardsrand(123) do - A = rand(elty,10,10) - A = A + A.' #symmetric! - b = rand(elty,10) - c = A \ b - b,A = LAPACK.sysv!('U',A,b) - @test b ≈ c - @test_throws DimensionMismatch LAPACK.sysv!('U',A,rand(elty,11)) - end + srand(123) + A = rand(elty,10,10) + A = A + A.' #symmetric! + b = rand(elty,10) + c = A \ b + b,A = LAPACK.sysv!('U',A,b) + @test b ≈ c + @test_throws DimensionMismatch LAPACK.sysv!('U',A,rand(elty,11)) end end @testset "hesv" begin @testset for elty in (ComplexF32, ComplexF64) - guardsrand(935) do - A = rand(elty,10,10) - A = A + A' #hermitian! - b = rand(elty,10) - c = A \ b - b,A = LAPACK.hesv!('U',A,b) - @test b ≈ c - @test_throws DimensionMismatch LAPACK.hesv!('U',A,rand(elty,11)) - A = rand(elty,10,10) - A = A + A' #hermitian! - b = rand(elty,10) - c = A \ b - b,A = LAPACK.hesv_rook!('U',A,b) - @test b ≈ c - @test_throws DimensionMismatch LAPACK.hesv_rook!('U',A,rand(elty,11)) - end + srand(935) + A = rand(elty,10,10) + A = A + A' #hermitian! + b = rand(elty,10) + c = A \ b + b,A = LAPACK.hesv!('U',A,b) + @test b ≈ c + @test_throws DimensionMismatch LAPACK.hesv!('U',A,rand(elty,11)) + A = rand(elty,10,10) + A = A + A' #hermitian! + b = rand(elty,10) + c = A \ b + b,A = LAPACK.hesv_rook!('U',A,b) + @test b ≈ c + @test_throws DimensionMismatch LAPACK.hesv_rook!('U',A,rand(elty,11)) end end diff --git a/test/linalg/lu.jl b/test/linalg/lu.jl index da0949491a32d..a19a72ed4b85e 100644 --- a/test/linalg/lu.jl +++ b/test/linalg/lu.jl @@ -195,6 +195,7 @@ dimg = randn(n)/2 end @testset "conversion" begin + srand(3) a = Tridiagonal(rand(9),rand(10),rand(9)) fa = Array(a) falu = lufact(fa) diff --git a/test/linalg/tridiag.jl b/test/linalg/tridiag.jl index bf87e97d309e0..8ba310548b2c6 100644 --- a/test/linalg/tridiag.jl +++ b/test/linalg/tridiag.jl @@ -14,318 +14,318 @@ function test_approx_eq_vecs(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error end end -guardsrand(123) do +@testset for elty in (Float32, Float64, ComplexF32, ComplexF64, Int) n = 12 #Size of matrix problem to test - @testset for elty in (Float32, Float64, ComplexF32, ComplexF64, Int) - if elty == Int - srand(61516384) - d = rand(1:100, n) - dl = -rand(0:10, n-1) - du = -rand(0:10, n-1) - v = rand(1:100, n) - B = rand(1:100, n, 2) - a = rand(1:100, n-1) - b = rand(1:100, n) - c = rand(1:100, n-1) - else - d = convert(Vector{elty}, 1 .+ randn(n)) - dl = convert(Vector{elty}, randn(n - 1)) - du = convert(Vector{elty}, randn(n - 1)) - v = convert(Vector{elty}, randn(n)) - B = convert(Matrix{elty}, randn(n, 2)) - a = convert(Vector{elty}, randn(n - 1)) - b = convert(Vector{elty}, randn(n)) - c = convert(Vector{elty}, randn(n - 1)) - if elty <: Complex - a += im*convert(Vector{elty}, randn(n - 1)) - b += im*convert(Vector{elty}, randn(n)) - c += im*convert(Vector{elty}, randn(n - 1)) - end + srand(123) + if elty == Int + srand(61516384) + d = rand(1:100, n) + dl = -rand(0:10, n-1) + du = -rand(0:10, n-1) + v = rand(1:100, n) + B = rand(1:100, n, 2) + a = rand(1:100, n-1) + b = rand(1:100, n) + c = rand(1:100, n-1) + else + d = convert(Vector{elty}, 1 .+ randn(n)) + dl = convert(Vector{elty}, randn(n - 1)) + du = convert(Vector{elty}, randn(n - 1)) + v = convert(Vector{elty}, randn(n)) + B = convert(Matrix{elty}, randn(n, 2)) + a = convert(Vector{elty}, randn(n - 1)) + b = convert(Vector{elty}, randn(n)) + c = convert(Vector{elty}, randn(n - 1)) + if elty <: Complex + a += im*convert(Vector{elty}, randn(n - 1)) + b += im*convert(Vector{elty}, randn(n)) + c += im*convert(Vector{elty}, randn(n - 1)) end - @test_throws DimensionMismatch SymTridiagonal(dl, ones(elty, n + 1)) - @test_throws ArgumentError SymTridiagonal(rand(n, n)) - @test_throws ArgumentError Tridiagonal(dl, dl, dl) - @test_throws ArgumentError convert(SymTridiagonal{elty}, Tridiagonal(dl, d, du)) + end + @test_throws DimensionMismatch SymTridiagonal(dl, ones(elty, n + 1)) + @test_throws ArgumentError SymTridiagonal(rand(n, n)) + @test_throws ArgumentError Tridiagonal(dl, dl, dl) + @test_throws ArgumentError convert(SymTridiagonal{elty}, Tridiagonal(dl, d, du)) - if elty != Int - @testset "issue #1490" begin - @test det(ones(elty,3,3)) ≈ zero(elty) atol=3*eps(real(one(elty))) - @test det(SymTridiagonal(elty[],elty[])) == one(elty) - end + if elty != Int + @testset "issue #1490" begin + @test det(ones(elty,3,3)) ≈ zero(elty) atol=3*eps(real(one(elty))) + @test det(SymTridiagonal(elty[],elty[])) == one(elty) end + end - @testset "constructor" begin - for (x, y) in ((d, dl), (GenericArray(d), GenericArray(dl))) - ST = (SymTridiagonal(x, y))::SymTridiagonal{elty, typeof(x)} - @test ST == Matrix(ST) - @test ST.dv === x - @test ST.ev === y - TT = (Tridiagonal(y, x, y))::Tridiagonal{elty, typeof(x)} - @test TT == Matrix(TT) - @test TT.dl === y - @test TT.d === x - @test TT.du === y - end - # enable when deprecations for 0.7 are dropped - # @test_throws MethodError SymTridiagonal(dv, GenericArray(ev)) - # @test_throws MethodError SymTridiagonal(GenericArray(dv), ev) - # @test_throws MethodError Tridiagonal(GenericArray(ev), dv, GenericArray(ev)) - # @test_throws MethodError Tridiagonal(ev, GenericArray(dv), ev) - end - @testset "interconversion of Tridiagonal and SymTridiagonal" begin - @test Tridiagonal(dl, d, dl) == SymTridiagonal(d, dl) - @test SymTridiagonal(d, dl) == Tridiagonal(dl, d, dl) - @test Tridiagonal(dl, d, du) + Tridiagonal(du, d, dl) == SymTridiagonal(2d, dl+du) - @test SymTridiagonal(d, dl) + Tridiagonal(dl, d, du) == Tridiagonal(dl + dl, d+d, dl+du) - @test convert(SymTridiagonal,Tridiagonal(SymTridiagonal(d, dl))) == SymTridiagonal(d, dl) - @test Array(convert(SymTridiagonal{ComplexF32},Tridiagonal(SymTridiagonal(d, dl)))) == convert(Matrix{ComplexF32}, SymTridiagonal(d, dl)) + @testset "constructor" begin + for (x, y) in ((d, dl), (GenericArray(d), GenericArray(dl))) + ST = (SymTridiagonal(x, y))::SymTridiagonal{elty, typeof(x)} + @test ST == Matrix(ST) + @test ST.dv === x + @test ST.ev === y + TT = (Tridiagonal(y, x, y))::Tridiagonal{elty, typeof(x)} + @test TT == Matrix(TT) + @test TT.dl === y + @test TT.d === x + @test TT.du === y end - @testset "tril/triu" begin - zerosd = fill!(similar(d), 0) - zerosdl = fill!(similar(dl), 0) - zerosdu = fill!(similar(du), 0) - @test_throws ArgumentError tril!(SymTridiagonal(d, dl), -n - 2) - @test_throws ArgumentError tril!(SymTridiagonal(d, dl), n) - @test_throws ArgumentError tril!(Tridiagonal(dl, d, du), -n - 2) - @test_throws ArgumentError tril!(Tridiagonal(dl, d, du), n) - @test tril(SymTridiagonal(d,dl)) == Tridiagonal(dl,d,zerosdl) - @test tril(SymTridiagonal(d,dl),1) == Tridiagonal(dl,d,dl) - @test tril(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,zerosd,zerosdl) - @test tril(SymTridiagonal(d,dl),-2) == Tridiagonal(zerosdl,zerosd,zerosdl) - @test tril(Tridiagonal(dl,d,du)) == Tridiagonal(dl,d,zerosdu) - @test tril(Tridiagonal(dl,d,du),1) == Tridiagonal(dl,d,du) - @test tril(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,zerosd,zerosdu) - @test tril(Tridiagonal(dl,d,du),-2) == Tridiagonal(zerosdl,zerosd,zerosdu) + # enable when deprecations for 0.7 are dropped + # @test_throws MethodError SymTridiagonal(dv, GenericArray(ev)) + # @test_throws MethodError SymTridiagonal(GenericArray(dv), ev) + # @test_throws MethodError Tridiagonal(GenericArray(ev), dv, GenericArray(ev)) + # @test_throws MethodError Tridiagonal(ev, GenericArray(dv), ev) + end + @testset "interconversion of Tridiagonal and SymTridiagonal" begin + @test Tridiagonal(dl, d, dl) == SymTridiagonal(d, dl) + @test SymTridiagonal(d, dl) == Tridiagonal(dl, d, dl) + @test Tridiagonal(dl, d, du) + Tridiagonal(du, d, dl) == SymTridiagonal(2d, dl+du) + @test SymTridiagonal(d, dl) + Tridiagonal(dl, d, du) == Tridiagonal(dl + dl, d+d, dl+du) + @test convert(SymTridiagonal,Tridiagonal(SymTridiagonal(d, dl))) == SymTridiagonal(d, dl) + @test Array(convert(SymTridiagonal{ComplexF32},Tridiagonal(SymTridiagonal(d, dl)))) == convert(Matrix{ComplexF32}, SymTridiagonal(d, dl)) + end + @testset "tril/triu" begin + zerosd = fill!(similar(d), 0) + zerosdl = fill!(similar(dl), 0) + zerosdu = fill!(similar(du), 0) + @test_throws ArgumentError tril!(SymTridiagonal(d, dl), -n - 2) + @test_throws ArgumentError tril!(SymTridiagonal(d, dl), n) + @test_throws ArgumentError tril!(Tridiagonal(dl, d, du), -n - 2) + @test_throws ArgumentError tril!(Tridiagonal(dl, d, du), n) + @test tril(SymTridiagonal(d,dl)) == Tridiagonal(dl,d,zerosdl) + @test tril(SymTridiagonal(d,dl),1) == Tridiagonal(dl,d,dl) + @test tril(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,zerosd,zerosdl) + @test tril(SymTridiagonal(d,dl),-2) == Tridiagonal(zerosdl,zerosd,zerosdl) + @test tril(Tridiagonal(dl,d,du)) == Tridiagonal(dl,d,zerosdu) + @test tril(Tridiagonal(dl,d,du),1) == Tridiagonal(dl,d,du) + @test tril(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,zerosd,zerosdu) + @test tril(Tridiagonal(dl,d,du),-2) == Tridiagonal(zerosdl,zerosd,zerosdu) - @test_throws ArgumentError triu!(SymTridiagonal(d, dl), -n) - @test_throws ArgumentError triu!(SymTridiagonal(d, dl), n + 2) - @test_throws ArgumentError triu!(Tridiagonal(dl, d, du), -n) - @test_throws ArgumentError triu!(Tridiagonal(dl, d, du), n + 2) - @test triu(SymTridiagonal(d,dl)) == Tridiagonal(zerosdl,d,dl) - @test triu(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,d,dl) - @test triu(SymTridiagonal(d,dl),1) == Tridiagonal(zerosdl,zerosd,dl) - @test triu(SymTridiagonal(d,dl),2) == Tridiagonal(zerosdl,zerosd,zerosdl) - @test triu(Tridiagonal(dl,d,du)) == Tridiagonal(zerosdl,d,du) - @test triu(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,d,du) - @test triu(Tridiagonal(dl,d,du),1) == Tridiagonal(zerosdl,zerosd,du) - @test triu(Tridiagonal(dl,d,du),2) == Tridiagonal(zerosdl,zerosd,zerosdu) + @test_throws ArgumentError triu!(SymTridiagonal(d, dl), -n) + @test_throws ArgumentError triu!(SymTridiagonal(d, dl), n + 2) + @test_throws ArgumentError triu!(Tridiagonal(dl, d, du), -n) + @test_throws ArgumentError triu!(Tridiagonal(dl, d, du), n + 2) + @test triu(SymTridiagonal(d,dl)) == Tridiagonal(zerosdl,d,dl) + @test triu(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,d,dl) + @test triu(SymTridiagonal(d,dl),1) == Tridiagonal(zerosdl,zerosd,dl) + @test triu(SymTridiagonal(d,dl),2) == Tridiagonal(zerosdl,zerosd,zerosdl) + @test triu(Tridiagonal(dl,d,du)) == Tridiagonal(zerosdl,d,du) + @test triu(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,d,du) + @test triu(Tridiagonal(dl,d,du),1) == Tridiagonal(zerosdl,zerosd,du) + @test triu(Tridiagonal(dl,d,du),2) == Tridiagonal(zerosdl,zerosd,zerosdu) - @test !istril(SymTridiagonal(d,dl)) - @test !istriu(SymTridiagonal(d,dl)) - @test istriu(Tridiagonal(zerosdl,d,du)) - @test istril(Tridiagonal(dl,d,zerosdu)) - end + @test !istril(SymTridiagonal(d,dl)) + @test !istriu(SymTridiagonal(d,dl)) + @test istriu(Tridiagonal(zerosdl,d,du)) + @test istril(Tridiagonal(dl,d,zerosdu)) + end - @testset for mat_type in (Tridiagonal, SymTridiagonal) - A = mat_type == Tridiagonal ? mat_type(dl, d, du) : mat_type(d, dl) - fA = map(elty <: Complex ? ComplexF64 : Float64, Array(A)) - @testset "similar, size, and copyto!" begin - B = similar(A) - @test size(B) == size(A) - if mat_type == Tridiagonal # doesn't work for SymTridiagonal yet - copyto!(B, A) - @test B == A - end - @test isa(similar(A), mat_type{elty}) - @test isa(similar(A, Int), mat_type{Int}) - @test isa(similar(A, (3, 2)), SparseMatrixCSC) - @test isa(similar(A, Int, (3, 2)), SparseMatrixCSC{Int}) - @test size(A, 3) == 1 - @test size(A, 1) == n - @test size(A) == (n, n) - @test_throws ArgumentError size(A, 0) - end - @testset "getindex" begin - @test_throws BoundsError A[n + 1, 1] - @test_throws BoundsError A[1, n + 1] - @test A[1, n] == convert(elty, 0.0) - @test A[1, 1] == d[1] + @testset for mat_type in (Tridiagonal, SymTridiagonal) + A = mat_type == Tridiagonal ? mat_type(dl, d, du) : mat_type(d, dl) + fA = map(elty <: Complex ? ComplexF64 : Float64, Array(A)) + @testset "similar, size, and copyto!" begin + B = similar(A) + @test size(B) == size(A) + if mat_type == Tridiagonal # doesn't work for SymTridiagonal yet + copyto!(B, A) + @test B == A end - @testset "setindex!" begin - @test_throws BoundsError A[n + 1, 1] = 0 # test bounds check - @test_throws BoundsError A[1, n + 1] = 0 # test bounds check - @test_throws ArgumentError A[1, 3] = 1 # test assignment off the main/sub/super diagonal - if mat_type == Tridiagonal - @test (A[3, 3] = A[3, 3]; A == fA) # test assignment on the main diagonal - @test (A[3, 2] = A[3, 2]; A == fA) # test assignment on the subdiagonal - @test (A[2, 3] = A[2, 3]; A == fA) # test assignment on the superdiagonal - @test ((A[1, 3] = 0) == 0; A == fA) # test zero assignment off the main/sub/super diagonal - else # mat_type is SymTridiagonal - @test ((A[3, 3] = A[3, 3]) == A[3, 3]; A == fA) # test assignment on the main diagonal - @test_throws ArgumentError A[3, 2] = 1 # test assignment on the subdiagonal - @test_throws ArgumentError A[2, 3] = 1 # test assignment on the superdiagonal - end + @test isa(similar(A), mat_type{elty}) + @test isa(similar(A, Int), mat_type{Int}) + @test isa(similar(A, (3, 2)), SparseMatrixCSC) + @test isa(similar(A, Int, (3, 2)), SparseMatrixCSC{Int}) + @test size(A, 3) == 1 + @test size(A, 1) == n + @test size(A) == (n, n) + @test_throws ArgumentError size(A, 0) + end + @testset "getindex" begin + @test_throws BoundsError A[n + 1, 1] + @test_throws BoundsError A[1, n + 1] + @test A[1, n] == convert(elty, 0.0) + @test A[1, 1] == d[1] + end + @testset "setindex!" begin + @test_throws BoundsError A[n + 1, 1] = 0 # test bounds check + @test_throws BoundsError A[1, n + 1] = 0 # test bounds check + @test_throws ArgumentError A[1, 3] = 1 # test assignment off the main/sub/super diagonal + if mat_type == Tridiagonal + @test (A[3, 3] = A[3, 3]; A == fA) # test assignment on the main diagonal + @test (A[3, 2] = A[3, 2]; A == fA) # test assignment on the subdiagonal + @test (A[2, 3] = A[2, 3]; A == fA) # test assignment on the superdiagonal + @test ((A[1, 3] = 0) == 0; A == fA) # test zero assignment off the main/sub/super diagonal + else # mat_type is SymTridiagonal + @test ((A[3, 3] = A[3, 3]) == A[3, 3]; A == fA) # test assignment on the main diagonal + @test_throws ArgumentError A[3, 2] = 1 # test assignment on the subdiagonal + @test_throws ArgumentError A[2, 3] = 1 # test assignment on the superdiagonal end - @testset "diag" begin - @test (@inferred diag(A))::typeof(d) == d - @test (@inferred diag(A, 0))::typeof(d) == d - @test (@inferred diag(A, 1))::typeof(d) == (mat_type == Tridiagonal ? du : dl) - @test (@inferred diag(A, -1))::typeof(d) == dl - @test (@inferred diag(A, n-1))::typeof(d) == zeros(elty, 1) - @test_throws ArgumentError diag(A, -n - 1) - @test_throws ArgumentError diag(A, n + 1) - GA = mat_type == Tridiagonal ? mat_type(GenericArray.((dl, d, du))...) : mat_type(GenericArray.((d, dl))...) - @test (@inferred diag(GA))::typeof(GenericArray(d)) == GenericArray(d) - @test (@inferred diag(GA, -1))::typeof(GenericArray(d)) == GenericArray(dl) + end + @testset "diag" begin + @test (@inferred diag(A))::typeof(d) == d + @test (@inferred diag(A, 0))::typeof(d) == d + @test (@inferred diag(A, 1))::typeof(d) == (mat_type == Tridiagonal ? du : dl) + @test (@inferred diag(A, -1))::typeof(d) == dl + @test (@inferred diag(A, n-1))::typeof(d) == zeros(elty, 1) + @test_throws ArgumentError diag(A, -n - 1) + @test_throws ArgumentError diag(A, n + 1) + GA = mat_type == Tridiagonal ? mat_type(GenericArray.((dl, d, du))...) : mat_type(GenericArray.((d, dl))...) + @test (@inferred diag(GA))::typeof(GenericArray(d)) == GenericArray(d) + @test (@inferred diag(GA, -1))::typeof(GenericArray(d)) == GenericArray(dl) + end + @testset "Idempotent tests" begin + for func in (conj, transpose, adjoint) + @test func(func(A)) == A end - @testset "Idempotent tests" begin - for func in (conj, transpose, adjoint) - @test func(func(A)) == A + end + if elty != Int + @testset "Simple unary functions" begin + for func in (det, inv) + @test func(A) ≈ func(fA) atol=n^2*sqrt(eps(real(one(elty)))) end end - if elty != Int - @testset "Simple unary functions" begin - for func in (det, inv) - @test func(A) ≈ func(fA) atol=n^2*sqrt(eps(real(one(elty)))) - end - end + end + ds = mat_type == Tridiagonal ? (dl, d, du) : (d, dl) + for f in (real, imag) + @test f(A)::mat_type == mat_type(map(f, ds)...) + end + if elty <: Real + for f in (round, trunc, floor, ceil) + fds = [f.(d) for d in ds] + @test f.(A)::mat_type == mat_type(fds...) + @test f.(Int, A)::mat_type == f.(Int, fA) end - ds = mat_type == Tridiagonal ? (dl, d, du) : (d, dl) - for f in (real, imag) - @test f(A)::mat_type == mat_type(map(f, ds)...) + end + fds = [abs.(d) for d in ds] + @test abs.(A)::mat_type == mat_type(fds...) + @testset "Multiplication with strided matrix/vector" begin + @test A*ones(n) ≈ Array(A)*ones(n) + @test A*ones(n, 2) ≈ Array(A)*ones(n, 2) + end + @testset "Binary operations" begin + B = mat_type == Tridiagonal ? mat_type(a, b, c) : mat_type(b, a) + fB = map(elty <: Complex ? ComplexF64 : Float64, Array(B)) + for op in (+, -, *) + @test Array(op(A, B)) ≈ op(fA, fB) end - if elty <: Real - for f in (round, trunc, floor, ceil) - fds = [f.(d) for d in ds] - @test f.(A)::mat_type == mat_type(fds...) - @test f.(Int, A)::mat_type == f.(Int, fA) - end + α = rand(elty) + @test Array(α*A) ≈ α*Array(A) + @test Array(A*α) ≈ Array(A)*α + @test Array(A/α) ≈ Array(A)/α + + @testset "Matmul with Triangular types" begin + @test A*Base.LinAlg.UnitUpperTriangular(Matrix(1.0I, n, n)) ≈ fA + @test A*Base.LinAlg.UnitLowerTriangular(Matrix(1.0I, n, n)) ≈ fA + @test A*UpperTriangular(Matrix(1.0I, n, n)) ≈ fA + @test A*LowerTriangular(Matrix(1.0I, n, n)) ≈ fA end - fds = [abs.(d) for d in ds] - @test abs.(A)::mat_type == mat_type(fds...) - @testset "Multiplication with strided matrix/vector" begin - @test A*ones(n) ≈ Array(A)*ones(n) - @test A*ones(n, 2) ≈ Array(A)*ones(n, 2) + @testset "mul! errors" begin + @test_throws DimensionMismatch Base.LinAlg.mul!(similar(fA),A,ones(elty,n,n+1)) + @test_throws DimensionMismatch Base.LinAlg.mul!(similar(fA),A,ones(elty,n+1,n)) + @test_throws DimensionMismatch Base.LinAlg.mul!(zeros(elty,n,n),B,ones(elty,n+1,n)) + @test_throws DimensionMismatch Base.LinAlg.mul!(zeros(elty,n+1,n),B,ones(elty,n,n)) + @test_throws DimensionMismatch Base.LinAlg.mul!(zeros(elty,n,n+1),B,ones(elty,n,n)) end - @testset "Binary operations" begin - B = mat_type == Tridiagonal ? mat_type(a, b, c) : mat_type(b, a) - fB = map(elty <: Complex ? ComplexF64 : Float64, Array(B)) - for op in (+, -, *) - @test Array(op(A, B)) ≈ op(fA, fB) - end - α = rand(elty) - @test Array(α*A) ≈ α*Array(A) - @test Array(A*α) ≈ Array(A)*α - @test Array(A/α) ≈ Array(A)/α - - @testset "Matmul with Triangular types" begin - @test A*Base.LinAlg.UnitUpperTriangular(Matrix(1.0I, n, n)) ≈ fA - @test A*Base.LinAlg.UnitLowerTriangular(Matrix(1.0I, n, n)) ≈ fA - @test A*UpperTriangular(Matrix(1.0I, n, n)) ≈ fA - @test A*LowerTriangular(Matrix(1.0I, n, n)) ≈ fA - end - @testset "mul! errors" begin - @test_throws DimensionMismatch Base.LinAlg.mul!(similar(fA),A,ones(elty,n,n+1)) - @test_throws DimensionMismatch Base.LinAlg.mul!(similar(fA),A,ones(elty,n+1,n)) - @test_throws DimensionMismatch Base.LinAlg.mul!(zeros(elty,n,n),B,ones(elty,n+1,n)) - @test_throws DimensionMismatch Base.LinAlg.mul!(zeros(elty,n+1,n),B,ones(elty,n,n)) - @test_throws DimensionMismatch Base.LinAlg.mul!(zeros(elty,n,n+1),B,ones(elty,n,n)) - end + end + if mat_type == SymTridiagonal + @testset "Tridiagonal/SymTridiagonal mixing ops" begin + B = convert(Tridiagonal{elty}, A) + @test B == A + @test B + A == A + B + @test B - A == A - B end - if mat_type == SymTridiagonal - @testset "Tridiagonal/SymTridiagonal mixing ops" begin - B = convert(Tridiagonal{elty}, A) - @test B == A - @test B + A == A + B - @test B - A == A - B - end - if elty <: Base.LinAlg.BlasReal - @testset "Eigensystems" begin - zero, infinity = convert(elty, 0), convert(elty, Inf) - @testset "stebz! and stein!" begin - w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) - evecs = LAPACK.stein!(b, a, w) + if elty <: Base.LinAlg.BlasReal + @testset "Eigensystems" begin + zero, infinity = convert(elty, 0), convert(elty, Inf) + @testset "stebz! and stein!" begin + w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) + evecs = LAPACK.stein!(b, a, w) - (e, v) = eig(SymTridiagonal(b, a)) - @test e ≈ w - test_approx_eq_vecs(v, evecs) - end - @testset "stein! call using iblock and isplit" begin - w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) - evecs = LAPACK.stein!(b, a, w, iblock, isplit) - test_approx_eq_vecs(v, evecs) - end - @testset "stegr! call with index range" begin - F = eigfact(SymTridiagonal(b, a),1:2) - fF = eigfact(Symmetric(Array(SymTridiagonal(b, a))),1:2) - Test.test_approx_eq_modphase(F[:vectors], fF[:vectors]) - @test F[:values] ≈ fF[:values] - end - @testset "stegr! call with value range" begin - F = eigfact(SymTridiagonal(b, a),0.0,1.0) - fF = eigfact(Symmetric(Array(SymTridiagonal(b, a))),0.0,1.0) - Test.test_approx_eq_modphase(F[:vectors], fF[:vectors]) - @test F[:values] ≈ fF[:values] - end - @testset "eigenvalues/eigenvectors of symmetric tridiagonal" begin - if elty === Float32 || elty === Float64 - DT, VT = @inferred eig(A) - @inferred eig(A, 2:4) - @inferred eig(A, 1.0, 2.0) - D, Vecs = eig(fA) - @test DT ≈ D - @test abs.(VT'Vecs) ≈ Matrix(elty(1)I, n, n) - Test.test_approx_eq_modphase(eigvecs(A), eigvecs(fA)) - #call to LAPACK.stein here - Test.test_approx_eq_modphase(eigvecs(A,eigvals(A)),eigvecs(A)) - elseif elty != Int - # check that undef is determined accurately even if type inference - # bails out due to the number of try/catch blocks in this code. - @test_throws UndefVarError fA - end + (e, v) = eig(SymTridiagonal(b, a)) + @test e ≈ w + test_approx_eq_vecs(v, evecs) + end + @testset "stein! call using iblock and isplit" begin + w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) + evecs = LAPACK.stein!(b, a, w, iblock, isplit) + test_approx_eq_vecs(v, evecs) + end + @testset "stegr! call with index range" begin + F = eigfact(SymTridiagonal(b, a),1:2) + fF = eigfact(Symmetric(Array(SymTridiagonal(b, a))),1:2) + Test.test_approx_eq_modphase(F[:vectors], fF[:vectors]) + @test F[:values] ≈ fF[:values] + end + @testset "stegr! call with value range" begin + F = eigfact(SymTridiagonal(b, a),0.0,1.0) + fF = eigfact(Symmetric(Array(SymTridiagonal(b, a))),0.0,1.0) + Test.test_approx_eq_modphase(F[:vectors], fF[:vectors]) + @test F[:values] ≈ fF[:values] + end + @testset "eigenvalues/eigenvectors of symmetric tridiagonal" begin + if elty === Float32 || elty === Float64 + DT, VT = @inferred eig(A) + @inferred eig(A, 2:4) + @inferred eig(A, 1.0, 2.0) + D, Vecs = eig(fA) + @test DT ≈ D + @test abs.(VT'Vecs) ≈ Matrix(elty(1)I, n, n) + Test.test_approx_eq_modphase(eigvecs(A), eigvecs(fA)) + #call to LAPACK.stein here + Test.test_approx_eq_modphase(eigvecs(A,eigvals(A)),eigvecs(A)) + elseif elty != Int + # check that undef is determined accurately even if type inference + # bails out due to the number of try/catch blocks in this code. + @test_throws UndefVarError fA end end end - if elty <: Real - Ts = SymTridiagonal(d, dl) - Fs = Array(Ts) - Tldlt = factorize(Ts) - @testset "symmetric tridiagonal" begin - @test_throws DimensionMismatch Tldlt\rand(elty,n+1) - @test size(Tldlt) == size(Ts) - if elty <: AbstractFloat - @test typeof(convert(Base.LinAlg.LDLt{Float32},Tldlt)) == - Base.LinAlg.LDLt{Float32,SymTridiagonal{elty,Vector{elty}}} - end - for vv in (copy(v), view(v, 1:n)) - invFsv = Fs\vv - x = Ts\vv - @test x ≈ invFsv - @test Array(Tldlt) ≈ Fs - end - - @testset "similar" begin - @test isa(similar(Ts), SymTridiagonal{elty}) - @test isa(similar(Ts, Int), SymTridiagonal{Int}) - @test isa(similar(Ts, (3, 2)), SparseMatrixCSC) - @test isa(similar(Ts, Int, (3, 2)), SparseMatrixCSC{Int}) - end + end + if elty <: Real + Ts = SymTridiagonal(d, dl) + Fs = Array(Ts) + Tldlt = factorize(Ts) + @testset "symmetric tridiagonal" begin + @test_throws DimensionMismatch Tldlt\rand(elty,n+1) + @test size(Tldlt) == size(Ts) + if elty <: AbstractFloat + @test typeof(convert(Base.LinAlg.LDLt{Float32},Tldlt)) == + Base.LinAlg.LDLt{Float32,SymTridiagonal{elty,Vector{elty}}} + end + for vv in (copy(v), view(v, 1:n)) + invFsv = Fs\vv + x = Ts\vv + @test x ≈ invFsv + @test Array(Tldlt) ≈ Fs + end - @test first(logabsdet(Tldlt)) ≈ first(logabsdet(Fs)) - @test last(logabsdet(Tldlt)) ≈ last(logabsdet(Fs)) - # just test that the det method exists. The numerical value of the - # determinant is unreliable - det(Tldlt) + @testset "similar" begin + @test isa(similar(Ts), SymTridiagonal{elty}) + @test isa(similar(Ts, Int), SymTridiagonal{Int}) + @test isa(similar(Ts, (3, 2)), SparseMatrixCSC) + @test isa(similar(Ts, Int, (3, 2)), SparseMatrixCSC{Int}) end + + @test first(logabsdet(Tldlt)) ≈ first(logabsdet(Fs)) + @test last(logabsdet(Tldlt)) ≈ last(logabsdet(Fs)) + # just test that the det method exists. The numerical value of the + # determinant is unreliable + det(Tldlt) end - else # mat_type is Tridiagonal - @testset "tridiagonal linear algebra" begin - for (BB, vv) in ((copy(B), copy(v)), (view(B, 1:n, 1), view(v, 1:n))) - @test A*vv ≈ fA*vv - invFv = fA\vv - @test A\vv ≈ invFv - # @test Base.solve(T,v) ≈ invFv - # @test Base.solve(T, B) ≈ F\B - Tlu = factorize(A) - x = Tlu\vv - @test x ≈ invFv - end + end + else # mat_type is Tridiagonal + @testset "tridiagonal linear algebra" begin + for (BB, vv) in ((copy(B), copy(v)), (view(B, 1:n, 1), view(v, 1:n))) + @test A*vv ≈ fA*vv + invFv = fA\vv + @test A\vv ≈ invFv + # @test Base.solve(T,v) ≈ invFv + # @test Base.solve(T, B) ≈ F\B + Tlu = factorize(A) + x = Tlu\vv + @test x ≈ invFv end end end end end + @testset "Issue 12068" begin @test SymTridiagonal([1, 2], [0])^3 == [1 0; 0 8] end diff --git a/test/random.jl b/test/random.jl index 9130c306de94e..75c78cbb43483 100644 --- a/test/random.jl +++ b/test/random.jl @@ -5,8 +5,8 @@ using Main.TestHelpers.OAs using Base.Random: Sampler, SamplerRangeFast, SamplerRangeInt -# Issue #6573 -guardsrand(0) do +@testset "Issue #6573" begin + srand(0) rand() x = rand(384) @test find(x .== rand()) == [] diff --git a/test/sorting.jl b/test/sorting.jl index 9626c5de9e9b1..002b3394b1372 100644 --- a/test/sorting.jl +++ b/test/sorting.jl @@ -226,73 +226,73 @@ function randn_with_nans(n,p) end @testset "advanced sorting" begin - guardsrand(0xdeadbeef) do - for n in [0:10; 100; 101; 1000; 1001] - local r - r = -5:5 - v = rand(r,n) - h = [sum(v .== x) for x in r] - - for rev in [false,true] - # insertion sort (stable) as reference - pi = sortperm(v, alg=InsertionSort, rev=rev) - @test pi == sortperm(float(v), alg=InsertionSort, rev=rev) - @test isperm(pi) - si = v[pi] - @test [sum(si .== x) for x in r] == h - @test issorted(si, rev=rev) - @test all(issorted,[pi[si.==x] for x in r]) - c = copy(v) - permute!(c, pi) - @test c == si - ipermute!(c, pi) - @test c == v - - # stable algorithms - for alg in [MergeSort] - p = sortperm(v, alg=alg, rev=rev) - @test p == sortperm(float(v), alg=alg, rev=rev) - @test p == pi - s = copy(v) - permute!(s, p) - @test s == si - ipermute!(s, p) - @test s == v - end - - # unstable algorithms - for alg in [QuickSort, PartialQuickSort(n)] - p = sortperm(v, alg=alg, rev=rev) - @test p == sortperm(float(v), alg=alg, rev=rev) - @test isperm(p) - @test v[p] == si - s = copy(v) - permute!(s, p) - @test s == si - ipermute!(s, p) - @test s == v - end + srand(0xdeadbeef) + for n in [0:10; 100; 101; 1000; 1001] + local r + r = -5:5 + v = rand(r,n) + h = [sum(v .== x) for x in r] + + for rev in [false,true] + # insertion sort (stable) as reference + pi = sortperm(v, alg=InsertionSort, rev=rev) + @test pi == sortperm(float(v), alg=InsertionSort, rev=rev) + @test isperm(pi) + si = v[pi] + @test [sum(si .== x) for x in r] == h + @test issorted(si, rev=rev) + @test all(issorted,[pi[si.==x] for x in r]) + c = copy(v) + permute!(c, pi) + @test c == si + ipermute!(c, pi) + @test c == v + + # stable algorithms + for alg in [MergeSort] + p = sortperm(v, alg=alg, rev=rev) + @test p == sortperm(float(v), alg=alg, rev=rev) + @test p == pi + s = copy(v) + permute!(s, p) + @test s == si + ipermute!(s, p) + @test s == v end - v = randn_with_nans(n,0.1) - # TODO: alg = PartialQuickSort(n) fails here - for alg in [InsertionSort, QuickSort, MergeSort], - rev in [false,true] - # test float sorting with NaNs - s = sort(v, alg=alg, rev=rev) - @test issorted(s, rev=rev) - @test reinterpret(UInt64,v[isnan.(v)]) == reinterpret(UInt64,s[isnan.(s)]) - - # test float permutation with NaNs + # unstable algorithms + for alg in [QuickSort, PartialQuickSort(n)] p = sortperm(v, alg=alg, rev=rev) + @test p == sortperm(float(v), alg=alg, rev=rev) @test isperm(p) - vp = v[p] - @test isequal(vp,s) - @test reinterpret(UInt64,vp) == reinterpret(UInt64,s) + @test v[p] == si + s = copy(v) + permute!(s, p) + @test s == si + ipermute!(s, p) + @test s == v end end + + v = randn_with_nans(n,0.1) + # TODO: alg = PartialQuickSort(n) fails here + for alg in [InsertionSort, QuickSort, MergeSort], + rev in [false,true] + # test float sorting with NaNs + s = sort(v, alg=alg, rev=rev) + @test issorted(s, rev=rev) + @test reinterpret(UInt64,v[isnan.(v)]) == reinterpret(UInt64,s[isnan.(s)]) + + # test float permutation with NaNs + p = sortperm(v, alg=alg, rev=rev) + @test isperm(p) + vp = v[p] + @test isequal(vp,s) + @test reinterpret(UInt64,vp) == reinterpret(UInt64,s) + end end end + @testset "sortperm" begin inds = [ 1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10, diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index 0ffab8462f96f..27445e9184d77 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -1691,21 +1691,20 @@ end end @testset "sparse matrix normestinv" begin - guardsrand(1234) do - Ac = sprandn(20,20,.5) + im* sprandn(20,20,.5) - Aci = ceil.(Int64, 100*sprand(20,20,.5)) + im*ceil.(Int64, sprand(20,20,.5)) - Ar = sprandn(20,20,.5) - Ari = ceil.(Int64, 100*Ar) - if Base.USE_GPL_LIBS - # NOTE: normestinv is probabilistic, so must be included in the guardsrand block - @test Base.SparseArrays.normestinv(Ac,3) ≈ norm(inv(Array(Ac)),1) atol=1e-4 - @test Base.SparseArrays.normestinv(Aci,3) ≈ norm(inv(Array(Aci)),1) atol=1e-4 - @test Base.SparseArrays.normestinv(Ar) ≈ norm(inv(Array(Ar)),1) atol=1e-4 - @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,0) - @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,21) - end - @test_throws DimensionMismatch Base.SparseArrays.normestinv(sprand(3,5,.9)) + srand(1234) + Ac = sprandn(20,20,.5) + im* sprandn(20,20,.5) + Aci = ceil.(Int64, 100*sprand(20,20,.5)) + im*ceil.(Int64, sprand(20,20,.5)) + Ar = sprandn(20,20,.5) + Ari = ceil.(Int64, 100*Ar) + if Base.USE_GPL_LIBS + # NOTE: normestinv is probabilistic, so requires a fixed seed (set above in srand(1234)) + @test Base.SparseArrays.normestinv(Ac,3) ≈ norm(inv(Array(Ac)),1) atol=1e-4 + @test Base.SparseArrays.normestinv(Aci,3) ≈ norm(inv(Array(Aci)),1) atol=1e-4 + @test Base.SparseArrays.normestinv(Ar) ≈ norm(inv(Array(Ar)),1) atol=1e-4 + @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,0) + @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,21) end + @test_throws DimensionMismatch Base.SparseArrays.normestinv(sprand(3,5,.9)) end @testset "issue #13008" begin @@ -1745,30 +1744,29 @@ end end @testset "factorization" begin - guardsrand(123) do - local A - A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2) - A = A + A' - @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A)))) - A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2) - A = A*A' - @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A)))) - A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) - A = A + A.' - @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A)))) - A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) - A = A*A.' - @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A)))) - @test factorize(triu(A)) == triu(A) - @test isa(factorize(triu(A)), UpperTriangular{Float64, SparseMatrixCSC{Float64, Int}}) - @test factorize(tril(A)) == tril(A) - @test isa(factorize(tril(A)), LowerTriangular{Float64, SparseMatrixCSC{Float64, Int}}) - @test !Base.USE_GPL_LIBS || factorize(A[:, 1:4])\ones(size(A, 1)) ≈ Array(A[:, 1:4])\ones(size(A, 1)) - @test_throws ErrorException chol(A) - @test_throws ErrorException lu(A) - @test_throws ErrorException eig(A) - @test_throws ErrorException inv(A) - end + srand(123) + local A + A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2) + A = A + A' + @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A)))) + A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2) + A = A*A' + @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A)))) + A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + A = A + A.' + @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A)))) + A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + A = A*A.' + @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A)))) + @test factorize(triu(A)) == triu(A) + @test isa(factorize(triu(A)), UpperTriangular{Float64, SparseMatrixCSC{Float64, Int}}) + @test factorize(tril(A)) == tril(A) + @test isa(factorize(tril(A)), LowerTriangular{Float64, SparseMatrixCSC{Float64, Int}}) + @test !Base.USE_GPL_LIBS || factorize(A[:, 1:4])\ones(size(A, 1)) ≈ Array(A[:, 1:4])\ones(size(A, 1)) + @test_throws ErrorException chol(A) + @test_throws ErrorException lu(A) + @test_throws ErrorException eig(A) + @test_throws ErrorException inv(A) end @testset "issue #13792, use sparse triangular solvers for sparse triangular solves" begin