Skip to content

Commit

Permalink
Merge pull request #133 from JuliaLinearAlgebra/an/stablenrg
Browse files Browse the repository at this point in the history 
Use StableRNGs for testing to avoid breakage on Julia nightly
  • Loading branch information
@andreasnoack
andreasnoack authored Jul 1, 2021
2 parents 12b92ab + de0e323 commit ab80c9e
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 47 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,9 @@ Arpack_jll = "~3.5"
julia = "1.3"

[extras]
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Random", "SparseArrays", "Test"]
test = ["SparseArrays", "StableRNGs", "Test"]
89 changes: 44 additions & 45 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Arpack
using Test, LinearAlgebra, SparseArrays, Random
using Test, LinearAlgebra, SparseArrays, StableRNGs

@testset "eigs" begin
Random.seed!(1234)
rng = StableRNG(1235)
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))
areal = sprandn(rng, n, n, 0.4)
breal = sprandn(rng, n, n, 0.4)
acmplx = complex.(sprandn(rng, n, n, 0.4), sprandn(rng, n, n, 0.4))
bcmplx = complex.(sprandn(rng, n, n, 0.4), sprandn(rng, n, n, 0.4))

testtol = 1e-6

Expand All @@ -28,17 +28,17 @@ using Test, LinearAlgebra, SparseArrays, Random

b = convert(SparseMatrixCSC{elty}, b)
bsym = copy(b') + b
bpd = b'*b
bpd = b'*b + I

(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, LinearAlgebra.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)
@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
Expand Down Expand Up @@ -81,7 +81,7 @@ using Test, LinearAlgebra, SparseArrays, Random
end

@testset "ArgumentErrors" begin
@test_throws ArgumentError eigs(rand(elty,2,2))
@test_throws ArgumentError eigs(rand(rng, elty, 2, 2))
@test_throws ArgumentError eigs(a, nev=-1)
@test_throws ArgumentError eigs(a, which=:Z)
@test_throws ArgumentError eigs(a, which=:BE)
Expand All @@ -90,32 +90,32 @@ using Test, LinearAlgebra, SparseArrays, Random
if elty == Float64
@test_throws ArgumentError eigs(a + copy(transpose(a)), which=:SI)
@test_throws ArgumentError eigs(a + copy(transpose(a)), which=:LI)
@test_throws ArgumentError eigs(a, sigma = rand(ComplexF32))
@test_throws ArgumentError eigs(a, sigma = rand(rng, ComplexF32))
end
end
end

@testset "Symmetric generalized with singular B" begin
Random.seed!(127)
rng = StableRNG(127)
n = 10
k = 3
A = randn(n,n); A = A'A
B = randn(n,k); B = B*B'
A = randn(rng, n, n); A = A'A
B = randn(rng, n, k); B = B*B'
@test sort(eigs(A, B, nev = k, sigma = 1.0, explicittransform=:none)[1]) ≈ sort(eigvals(A, B); by=abs)[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
]
@testset "Problematic example from #6965A" begin
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
]
d, = eigs(A6965,which=:SM,nev=2,ncv=4,tol=eps())
@test d[1] ≈ 2.5346936860350002
@test real(d[2]) ≈ 2.6159972444834976
Expand Down Expand Up @@ -148,18 +148,17 @@ function LinearAlgebra.mul!(rho2::StridedVector{T},Phi::CPM{T},rho::StridedVecto
return copyto!(rho2,rho1)
end

let
# Generate random isometry
@testset "Test random isometry" begin
(Q, R) = qr(randn(100, 50))
Q = reshape(Array(Q), (50, 2, 50))
# Construct trace-preserving completely positive map from this
Phi = CPM(copy(Q))
(d,v,nconv,numiter,numop,resid) = eigs(Phi,nev=1,which=:LM)
(d,v,nconv,numiter,numop,resid) = eigs(Phi, nev=1, which=:LM)
# Properties: largest eigenvalue should be 1, largest eigenvector, when reshaped as matrix
# should be a Hermitian positive definite matrix (up to an arbitrary phase)

@test d[1] ≈ 1. # largest eigenvalue should be 1.
v = reshape(v,(50,50)) # reshape to matrix
v = reshape(v, (50, 50)) # reshape to matrix
v /= tr(v) # factor out arbitrary phase
@test norm(imag(v)) ≈ 0. # it should be real
v = real(v)
Expand All @@ -169,8 +168,8 @@ let
@test isposdef(v)

# Repeat with starting vector
(d2,v2,nconv2,numiter2,numop2,resid2) = eigs(Phi,nev=1,which=:LM,v0=reshape(v,(2500,)))
v2 = reshape(v2,(50,50))
(d2, v2, nconv2, numiter2, numop2, resid2) = eigs(Phi, nev=1, which=:LM, v0=reshape(v, (2500,)))
v2 = reshape(v2, (50,50))
v2 /= tr(v2)
@test numiter2 < numiter
@test v ≈ v2
Expand Down Expand Up @@ -210,16 +209,16 @@ end
@test S4[1].S ≈ [34.0, 6.0]
end
@testset "passing guess for Krylov vectors" begin
S1 = svds(A, nsv = 2, v0=rand(eltype(A),size(A,2)))
S1 = svds(A, nsv = 2, v0=rand(eltype(A), size(A,2)))
@test S1[1].S ≈ S2.S[1:2]
end

@test_throws ArgumentError svds(A,nsv=0)
@test_throws ArgumentError svds(A,nsv=20)
@test_throws DimensionMismatch svds(A,nsv=2,v0=rand(size(A,2)+1))
@test_throws ArgumentError svds(A, nsv=0)
@test_throws ArgumentError svds(A, nsv=20)
@test_throws DimensionMismatch svds(A, nsv=2, v0=rand(size(A,2) + 1))

@testset "Orthogonal vectors with repeated singular values $i times. Issue 16608" for i in 2:3
rng = MersenneTwister(126) # Fragile to compute repeated values without blocking so we set the seed
rng = StableRNG(126) # Fragile to compute repeated values without blocking so we set the seed
v0 = randn(rng, 20)
d = sort(rand(rng, 20), rev = true)
for j in 2:i
Expand Down Expand Up @@ -262,7 +261,7 @@ end
@test s1_right ≈ s2_right
end
@testset "passing guess for Krylov vectors" begin
S1 = svds(A, nsv = 2, v0=rand(eltype(A),size(A,2)))
S1 = svds(A, nsv = 2, v0=rand(eltype(A), size(A,2)))
@test S1[1].S ≈ S2.S[1:2]
end

Expand Down Expand Up @@ -295,16 +294,16 @@ LinearAlgebra.adjoint(A::MyOp) = MyOp(adjoint(A.mat))
end

@testset "low rank" begin
Random.seed!(123)
rng = StableRNG(123)
@testset "$T coefficients" for T in [Float64, Complex{Float64}]
@testset "rank $r" for r in [2, 5, 10]
m, n = 3*r, 4*r
nsv = 2*r

FU = qr(randn(T, m, r))
FU = qr(randn(rng, T, m, r))
U = Matrix(FU.Q)
S = 0.1 .+ sort(rand(r), rev=true)
FV = qr(randn(T, n, r))
S = 0.1 .+ sort(rand(rng, r), rev=true)
FV = qr(randn(rng, T, n, r))
V = Matrix(FV.Q)

A = U*Diagonal(S)*V'
Expand All @@ -325,9 +324,9 @@ end
end


# Problematic examples from #41
@test all(Matrix(svds([1. 0.; 0. 0.],nsv=1)[1]) ≈ [1. 0.; 0. 0.] for i in 1:10)
let A = [1. 0. 0.; 0. 0. 0.; 0. 0. 0.]
@testset "Problematic examples from #41" begin
@test all(Matrix(svds([1. 0.; 0. 0.],nsv=1)[1]) ≈ [1. 0.; 0. 0.] for i in 1:10)
A = [1. 0. 0.; 0. 0. 0.; 0. 0. 0.]
U,s,V = svds(A,nsv=2)[1]
@test U*Diagonal(s)*V' ≈ A atol=1e-7
@test U'U ≈ I
Expand Down

0 comments on commit ab80c9e

Please sign in to comment.