Merge pull request #218 from nep-pack/gallery_matrices

Gallery matrices closes #188

src/gallery_extra/basic_random_examples.jl

@@ -1,8 +1,8 @@
 # A delay eigenvalue problem
 function dep0(n::Int=5)
-    Random.seed!(0) # reset the random seed
-    A0 = randn(n,n)
-    A1 = randn(n,n)
+    msws_rng = MSWS_RNG()
+    A0 = gen_rng_mat(msws_rng,n,n)
+    A1 = gen_rng_mat(msws_rng,n,n)
     tau = 1.0
     nep = DEP([A0,A1],[0,tau])
     return nep
@@ -11,9 +11,9 @@ end
 
 # A delay eigenvalue problem with sparse matrices
 function dep0_sparse(n::Integer=100, p::Real=0.25)
-    Random.seed!(0) # reset the random seed
-    A0 = sparse(1:n,1:n,rand(n))+sprand(n,n,p)
-    A1 = sparse(1:n,1:n,rand(n))+sprand(n,n,p)
+    msws_rng = MSWS_RNG()
+    A0 = sparse(1:n,1:n,vec(gen_rng_mat(msws_rng,n,1)))+gen_rng_spmat(msws_rng,n,n,p)
+    A1 = sparse(1:n,1:n,vec(gen_rng_mat(msws_rng,n,1)))+gen_rng_spmat(msws_rng,n,n,p)
     tau = 1.0
     nep = DEP([A0,A1],[0,tau])
     return nep
@@ -22,10 +22,10 @@ end
 
 # A delay eigenvalue problem with sparse tridiagonal matrices
 function dep0_tridiag(n::Integer=100)
-    Random.seed!(1) # reset the random seed
+    msws_rng = MSWS_RNG()
     K = [1:n;2:n;1:n-1]; J=[1:n;1:n-1;2:n] # sparsity pattern of tridiag matrix
-    A0 = sparse(K, J, rand(3*n-2))
-    A1 = sparse(K, J, rand(3*n-2))
+    A0 = sparse(K, J, vec(gen_rng_mat(msws_rng,3*n-2,1)))
+    A1 = sparse(K, J, vec(gen_rng_mat(msws_rng,3*n-2,1)))
     tau = 1.0
     nep = DEP([A0,A1],[0,tau])
     return nep
@@ -34,10 +34,10 @@ end
 
 # A polynomial eigenvalue problem
 function pep0(n::Integer=200)
-    Random.seed!(0)
-    A0=randn(n,n)
-    A1=randn(n,n)
-    A2=randn(n,n)
+    msws_rng = MSWS_RNG()
+    A0=gen_rng_mat(msws_rng,n,n)
+    A1=gen_rng_mat(msws_rng,n,n)
+    A2=gen_rng_mat(msws_rng,n,n)
     A=[A0,A1,A2]
     nep=PEP(A)
     return nep
@@ -46,23 +46,81 @@ end
 
 # A polynomial eigenvalue problem with symmetric matrices
 function pep0_sym(n::Integer=200)
-    Random.seed!(0)
-    A0 = Symmetric(randn(n,n))
-    A1 = Symmetric(randn(n,n))
-    A2 = Symmetric(randn(n,n))
-    A = [A0]#, A1, A2]
+    msws_rng = MSWS_RNG()
+    A0 = Symmetric(gen_rng_mat(msws_rng,n,n))
+    A1 = Symmetric(gen_rng_mat(msws_rng,n,n))
+    A2 = Symmetric(gen_rng_mat(msws_rng,n,n))
+    A = [A0, A1, A2]
     nep = PEP(A)
     return nep
 end
 
 
 # A sparse polynomial eigenvalue problem
 function pep0_sparse(n::Integer=200, p::Real=0.03)
-    Random.seed!(0)
-    A0=sprandn(n,n,p)
-    A1=sprandn(n,n,p)
-    A2=sprandn(n,n,p)
+    msws_rng = MSWS_RNG()
+    A0=gen_rng_spmat(msws_rng,n,n,p)
+    A1=gen_rng_spmat(msws_rng,n,n,p)
+    A2=gen_rng_spmat(msws_rng,n,n,p)
     A=[A0,A1,A2]
     nep=PEP(A)
     return nep
 end
+
+
+# Helper for random matrices, stable over releases
+# B. Widynski, Middle Square Weyl Sequence RNG, arXiv 1704.00358
+mutable struct MSWS_RNG
+    x::UInt128
+    w::UInt128
+    s::UInt128
+    function MSWS_RNG(seed::UInt128 = zero(UInt128))
+        base = 0x9ef09a97ac0f9ecaef01c4f2db0958c9
+        s = (seed << 1) + base
+        x = 0x1de568e1a1ca1b593cbf13f7407cf43e
+        w = 0xd4ac5c288559e14a5fafc1b7df9f9e0e
+        return new(x, w, s)
+    end
+end
+
+function gen_rng_int(rng::MSWS_RNG)
+    rng.x *= rng.x
+    rng.x += (rng.w += rng.s)
+    rng.x = (rng.x>>64) | (rng.x<<64)
+    return UInt64(rng.x & typemax(UInt64))
+end
+
+function gen_rng_float(rng::MSWS_RNG)
+    return Float64(gen_rng_int(rng)/typemax(UInt64))
+end
+
+function gen_rng_mat(rng::MSWS_RNG, n::Int64, m::Int64)
+    A = zeros(Float64,n,m)
+    for c = 1:m
+        for r = 1:n
+            A[r,c] = 1-2*gen_rng_float(rng)
+        end
+    end
+    return A
+end
+
+function gen_rng_spmat(rng::MSWS_RNG, n::Int64, m::Int64, p::Real)
+    nonzeros = round(p*m*n)
+    dict = Dict{Tuple{Int64,Int64},Float64}()
+    for i = 1:nonzeros
+        r = mod(gen_rng_int(rng),n) + 1
+        c = mod(gen_rng_int(rng),m) + 1
+        dict[r,c] = 1-2*gen_rng_float(rng)
+    end
+    idxes = collect(keys(dict))
+    nnzs = length(idxes)
+    r = zeros(Int64,nnzs)
+    c = zeros(Int64,nnzs)
+    vals = zeros(Float64,nnzs)
+    for i = 1:nnzs
+        r[i] = idxes[i][1]
+        c[i] = idxes[i][2]
+        vals[i] = dict[idxes[i]]
+    end
+    return sparse(r,c,vals,n,m)
+end

src/gallery_extra/gallery_examples.jl

@@ -106,8 +106,8 @@ end
 # A quadratic eigenvalue problem with chosen eigenvalues
 function qep_fixed_eig(n::Integer=5, E::AbstractVecOrMat=NaN*ones(2*n))
     if any(isnan.(E))
-        Random.seed!(0) # reset the random seed
-        E=randn(2*n)
+        msws_rng = MSWS_RNG()
+        E=vec(gen_rng_mat(msws_rng,2*n,1))
     end
     A1 = diagm(0 => E[1:n])
     A2 = diagm(0 => E[n+1:2*n])

src/gallery_extra/periodic_dde.jl

@@ -315,11 +315,11 @@ function periodic_dde_gallery(::Type{PeriodicDDE_NEP}; name::String="mathieu",n=
         # 5.73989+0.732386im
         # 5.73989-0.732386im
 
-        Random.seed!(0);
-        A0=sprandn(n,n,0.3)-one(TT)*I
-        A1=sprandn(n,n,0.3)-one(TT)*I
-        B0=sprandn(n,n,0.3)-one(TT)*I
-        B1=sprandn(n,n,0.3)-one(TT)*I
+        msws_rng = MSWS_RNG()
+        A0=gen_rng_spmat(msws_rng,n,n,0.3)-one(TT)*I
+        A1=gen_rng_spmat(msws_rng,n,n,0.3)-one(TT)*I
+        B0=gen_rng_spmat(msws_rng,n,n,0.3)-one(TT)*I
+        B1=gen_rng_spmat(msws_rng,n,n,0.3)-one(TT)*I
         tau=2;
         A=t-> A0+cos(pi*t)*A1;
         B=t-> B0+exp(0.01*sin(pi*t))*B1;

test/beyn.jl

@@ -3,16 +3,19 @@
 using NonlinearEigenproblemsTest
 using NonlinearEigenproblems
 using Test
+using Random
 using LinearAlgebra
 
 
 @testset "Beyn contour" begin
+
     nep=nep_gallery("dep0")
     @bench @testset "disk at origin" begin
+        Random.seed!(0);
 
-        λ,v=contour_beyn(nep,logger=displaylevel,radius=1,neigs=2,sanity_check=false)
+        λ,v=contour_beyn(nep,logger=displaylevel,radius=1,neigs=1,sanity_check=false)
 
-        for i = 1:2
+        for i = 1:size(λ,1)
             @info "$i: $(λ[i])"
             M=compute_Mder(nep,λ[i])
             minimum(svdvals(M))
@@ -21,19 +24,19 @@ using LinearAlgebra
         end
 
 
-        λ,v=contour_beyn(nep,logger=displaylevel,radius=1,k=2,N=100,sanity_check=false)
+        λ,v=contour_beyn(nep,logger=displaylevel,radius=0.4,k=2,N=100,sanity_check=false)
         M=compute_Mder(nep,λ[1])
         minimum(svdvals(M))
         @test minimum(svdvals(M))<eps()*1000
 
     end
     @bench @testset "shifted disk" begin
 
-        λ,v=contour_beyn(nep,logger=displaylevel,σ=-0.2,radius=1.5,
+        λ,v=contour_beyn(nep,logger=displaylevel,σ=0.2,radius=1.0,
                          neigs=4,sanity_check=false)
 
-        @test size(λ,1)==4
-        for i = 1:4
+        @test size(λ,1)==3
+        for i = 1:3
             @info "$i: $(λ[i])"
             M=compute_Mder(nep,λ[i])
             minimum(svdvals(M))