use a better seed than 0

Author: stevengj

src/IterativeSolvers.jl

@@ -6,6 +6,10 @@ eigensystems, and singular value problems.
 """
 module IterativeSolvers
 
+# deterministic seed for pseudo-random numbers
+# (taken from the IterativeSolvers.jl UUID)
+const seed = [0x42fd0dbc, 0xa9815370, 0x80f2aaf5, 0x04508153]
+
 include("common.jl")
 include("orthogonalize.jl")
 include("history.jl")

src/bicgstabl.jl

@@ -29,7 +29,7 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2;
                              max_mv_products = size(A, 2),
                              abstol::Real = zero(real(eltype(b))),
                              reltol::Real = sqrt(eps(real(eltype(b)))),
-                             rng::AbstractRNG=MersenneTwister(0),
+                             rng::AbstractRNG=MersenneTwister(seed),
                              initial_zero = false)
     T = eltype(x)
     n = size(A, 1)

src/idrs.jl

@@ -82,7 +82,7 @@ end
 
 function idrs_method!(log::ConvergenceHistory, X, A, C::T,
     s::Number, Pl::precT, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false,
-    rng::AbstractRNG=MersenneTwister(0)) where {T, precT}
+    rng::AbstractRNG=MersenneTwister(seed)) where {T, precT}
 
     verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm")
     R = C - A*X

src/lobpcg.jl

@@ -787,7 +787,7 @@ Finds the `nev` extremal eigenvalues and their corresponding eigenvectors satisf
 function lobpcg(A, largest::Bool, nev::Int; kwargs...)
     lobpcg(A, nothing, largest, nev; kwargs...)
 end
-function lobpcg(A, B, largest::Bool, nev::Int; rng::AbstractRNG=MersenneTwister(0), kwargs...)
+function lobpcg(A, B, largest::Bool, nev::Int; rng::AbstractRNG=MersenneTwister(seed), kwargs...)
     lobpcg(A, B, largest, rand(rng, eltype(A), size(A, 1), nev); not_zeros=true, rng=rng, kwargs...)
 end
 
@@ -827,7 +827,7 @@ function lobpcg(A, largest::Bool, X0; kwargs...)
 end
 function lobpcg(A, B, largest, X0;
                 not_zeros=false, log=false, P=nothing, maxiter=200,
-                rng::AbstractRNG=MersenneTwister(0),
+                rng::AbstractRNG=MersenneTwister(seed),
                 C=nothing, tol::Real=default_tolerance(eltype(X0)))
     X = copy(X0)
     n = size(X, 1)
@@ -866,7 +866,7 @@ end
 
 """
 function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=false,
-                 rng::AbstractRNG=MersenneTwister(0),
+                 rng::AbstractRNG=MersenneTwister(seed),
                  tol::Real=default_tolerance(eltype(iterator.XBlocks.block)))
     X = iterator.XBlocks.block
     iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
@@ -931,7 +931,7 @@ function lobpcg(A, largest::Bool, X0, nev::Int; kwargs...)
 end
 function lobpcg(A, B, largest::Bool, X0, nev::Int;
                 not_zeros=false, log=false, P=nothing, maxiter=200,
-                rng::AbstractRNG=MersenneTwister(0),
+                rng::AbstractRNG=MersenneTwister(seed),
                 C=nothing, tol::Real=default_tolerance(eltype(X0)))
     n = size(X0, 1)
     sizeX = size(X0, 2)

src/simple.jl

@@ -62,7 +62,7 @@ See [`powm!`](@ref). Calls `powm!(B, x0; kwargs...)` with
 `x0` initialized as a random, complex unit vector using
 `rng::AbstractRNG`.
 """
-function powm(B; rng::AbstractRNG=MersenneTwister(0), kwargs...)
+function powm(B; rng::AbstractRNG=MersenneTwister(seed), kwargs...)
     x0 = rand(rng, Complex{real(eltype(B))}, size(B, 1))
     rmul!(x0, one(eltype(B)) / norm(x0))
     powm!(B, x0; kwargs...)
@@ -169,7 +169,7 @@ Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), 50, ismutating = true)
 λ, x = invpowm(Fmap, shift = σ, tol = 1e-4, maxiter = 200)
 ```
 """
-function invpowm(B; rng::AbstractRNG=MersenneTwister(0), kwargs...)
+function invpowm(B; rng::AbstractRNG=MersenneTwister(seed), kwargs...)
     x0 = rand(rng, Complex{real(eltype(B))}, size(B, 1))
     rmul!(x0, one(eltype(B)) / norm(x0))
     invpowm!(B, x0; kwargs...)

src/svdl.jl

@@ -176,7 +176,7 @@ end
 #########################
 
 function svdl_method!(log::ConvergenceHistory, A, l::Int=min(6, size(A,1)); k::Int=2l,
-    j::Int=l, rng::AbstractRNG=MersenneTwister(l), v0::AbstractVector = Vector{eltype(A)}(randn(rng, size(A, 2))) |> x->rmul!(x, inv(norm(x))),
+    j::Int=l, rng::AbstractRNG=MersenneTwister(seed), v0::AbstractVector = Vector{eltype(A)}(randn(rng, size(A, 2))) |> x->rmul!(x, inv(norm(x))),
     maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps(),
     verbose::Bool=false, method::Symbol=:ritz, vecs=:none, dolock::Bool=false)