use deterministic RNG initial state

Project.toml

@@ -1,6 +1,6 @@
 name = "IterativeSolvers"
 uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
-version = "0.9.2"
+version = "0.9.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -11,4 +11,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 RecipesBase = "0.6, 0.7, 0.8, 1.0"
-julia = "1.3"
+julia = "1.6"

docs/src/eigenproblems/lobpcg.md

@@ -13,6 +13,10 @@ lobpcg!
 
 A `LOBPCGIterator` is created to pre-allocate all the memory required by the method using the constructor `LOBPCGIterator(A, B, largest, X, P, C)` where `A` and `B` are the matrices from the generalized eigenvalue problem, `largest` indicates if the problem is a maximum or minimum eigenvalue problem, `X` is the initial eigenbasis, randomly sampled if not input, where `size(X, 2)` is the block size `bs`. `P` is the preconditioner, `nothing` by default, and `C` is the constraints matrix. The desired `k` eigenvalues are found `bs` at a time.
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument.
 
 ## References
 Implementation is based on [^Knyazev1993] and [^Scipy].

docs/src/linear_systems/bicgstabl.md

@@ -15,11 +15,16 @@ The method is based on the original article [^Sleijpen1993], but does not implem
 
 The `r` and `u` factors are pre-allocated as matrices of size $n \times (l + 1)$, so that BLAS2 methods can be used. Also the random shadow residual is pre-allocated as a vector. Hence the storage costs are approximately $2l + 3$ vectors.
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument.
+
 !!! tip
  BiCGStabl(l) can be used as an [iterator](@ref Iterators).
 
-[^Sleijpen1993]: 
+[^Sleijpen1993]:
 
- Sleijpen, Gerard LG, and Diederik R. Fokkema. "BiCGstab(l) for 
- linear equations involving unsymmetric matrices with complex spectrum." 
+ Sleijpen, Gerard LG, and Diederik R. Fokkema. "BiCGstab(l) for
+ linear equations involving unsymmetric matrices with complex spectrum."
  Electronic Transactions on Numerical Analysis 1.11 (1993): 2000.

docs/src/linear_systems/idrs.md

@@ -19,3 +19,7 @@ and the IDR chapter in [^Meurant2020].
 [^VanGijzen2011]: Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits Bi-orthogonality Properties. M. B. van Gijzen and P. Sonneveld ACM Trans. Math. Software, Vol. 38, No. 1, pp. 5:1-5:19, 2011
 [^Meurant2020]: The IDR family. G. Meurant and J. Duintjer Tebbens. In: Krylov Methods for Nonsymmetric Linear Systems. Springer Series in Computational Mathematics, vol 57. Springer, 2020. [doi:10.1007/978-3-030-55251-0_10](https://doi.org/10.1007/978-3-030-55251-0_10)
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument.

docs/src/svd/svdl.md

@@ -14,4 +14,10 @@ The implementation of thick restarting follows closely that of SLEPc as describe
 
 The singular vectors are computed directly by forming the Ritz vectors from the product of the Lanczos vectors `L.P`/`L.Q` and the singular vectors of `L.B`. Additional accuracy in the singular triples can be obtained using inverse iteration.
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument, or by passing an initial `v0` vector
+directly.
+
 [^Hernandez2008]: Vicente Hernández, José E. Román, and Andrés Tomás. "A robust and efficient parallel SVD solver based on restarted Lanczos bidiagonalization." Electronic Transactions on Numerical Analysis 31 (2008): 68-85.

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

@@ -1,5 +1,5 @@
 export bicgstabl, bicgstabl!, bicgstabl_iterator, bicgstabl_iterator!, BiCGStabIterable
-using Printf
+using Printf, Random
 import Base: iterate
 
 mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, smallMatT <: AbstractMatrix, realT <: Real, scalarT <: Number}
@@ -29,13 +29,14 @@ 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(seed),
  initial_zero = false)
  T = eltype(x)
  n = size(A, 1)
  mv_products = 0
 
  # Large vectors.
- r_shadow = rand(T, n)
+ r_shadow = rand(rng, T, n)
  rs = Matrix{T}(undef, n, l + 1)
  us = zeros(T, n, l + 1)
 
@@ -156,6 +157,7 @@ bicgstabl(A, b, l = 2; kwargs...) = bicgstabl!(zerox(A, b), A, b, l; initial_zer
 - `max_mv_products::Int = size(A, 2)`: maximum number of matrix vector products.
 For BiCGStab(l) this is a less dubious term than "number of iterations";
 - `Pl = Identity()`: left preconditioner of the method;
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 - `abstol::Real = zero(real(eltype(b)))`,
  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
  tolerance for the stopping condition

src/idrs.jl

@@ -32,6 +32,7 @@ shadow space.
  is the residual in the `k`th iteration;
 - `maxiter::Int = size(A, 2)`: maximum number of iterations;
 - `log::Bool`: keep track of the residual norm in each iteration;
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 - `verbose::Bool`: print convergence information during the iterations.
 
 # Return values
@@ -80,8 +81,8 @@ end
 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
- ) where {T, precT}
+ s::Number, Pl::precT, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false,
+ rng::AbstractRNG=MersenneTwister(seed)) where {T, precT}
 
  verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm")
  R = C - A*X
@@ -102,7 +103,7 @@ function idrs_method!(log::ConvergenceHistory, X, A, C::T,
 
  Z = zero(C)
 
- P = T[rand!(copy(C)) for k in 1:s]
+ P = T[rand!(rng, copy(C)) for k in 1:s]
  U = T[copy(Z) for k in 1:s]
  G = T[copy(Z) for k in 1:s]
  Q = copy(Z)

src/lobpcg.jl

@@ -787,8 +787,8 @@ 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; kwargs...)
- lobpcg(A, B, largest, rand(eltype(A), size(A, 1), nev); not_zeros=true, 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
 
 """
@@ -804,6 +804,7 @@ end
 ## Keywords
 
 - `not_zeros`: default is `false`. If `true`, `X0` will be assumed to not have any all-zeros column.
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 
 - `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
  states; if `false` only `results.trace` will be empty;
@@ -826,6 +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(seed),
  C=nothing, tol::Real=default_tolerance(eltype(X0)))
  X = copy(X0)
  n = size(X, 1)
@@ -835,7 +837,7 @@ function lobpcg(A, B, largest, X0;
 
  iterator = LOBPCGIterator(A, B, largest, X, P, C)
 
- return lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
+ return lobpcg!(iterator, rng=rng, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
 end
 
 """
@@ -849,6 +851,7 @@ end
 ## Keywords
 
 - `not_zeros`: default is `false`. If `true`, the initial Ritz vectors will be assumed to not have any all-zeros column.
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 
 - `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
  states; if `false` only `results.trace` will be empty;
@@ -863,21 +866,22 @@ end
 
 """
 function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=false,
+ rng::AbstractRNG=MersenneTwister(seed),
  tol::Real=default_tolerance(eltype(iterator.XBlocks.block)))
  X = iterator.XBlocks.block
  iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
  if !not_zeros
  @views for j in 1:size(X,2)
  if all(x -> x==0, X[:, j])
- @inbounds X[:,j] .= rand.()
+ @inbounds X[:,j] .= rand.(rng)
  end
  end
  iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
  end
  n = size(X, 1)
  sizeX = size(X, 2)
  iterator.iteration[] = 1
- while iterator.iteration[] <= maxiter 
+ while iterator.iteration[] <= maxiter
  state = iterator(tol, log)
  if log
  push!(iterator.trace, state)
@@ -927,6 +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(seed),
  C=nothing, tol::Real=default_tolerance(eltype(X0)))
  n = size(X0, 1)
  sizeX = size(X0, 2)
@@ -938,22 +943,22 @@ function lobpcg(A, B, largest::Bool, X0, nev::Int;
  iterator = LOBPCGIterator(A, B, largest, X, nev, P, C)
 
  r = EmptyLOBPCGResults(X, nev, tol, maxiter)
- rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
+ rnext = lobpcg!(iterator, rng=rng, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
  append!(r, rnext, 0)
  converged_x = sizeX
  while converged_x < nev
  @views if nev-converged_x < sizeX
  cutoff = sizeX-(nev-converged_x)
  update!(iterator.constr!, iterator.XBlocks.block[:, 1:cutoff], iterator.XBlocks.B_block[:, 1:cutoff])
  X[:, 1:sizeX-cutoff] .= X[:, cutoff+1:sizeX]
- rand!(X[:, cutoff+1:sizeX])
- rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
+ rand!(rng, X[:, cutoff+1:sizeX])
+ rnext = lobpcg!(iterator, rng=rng, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
  append!(r, rnext, converged_x, sizeX-cutoff)
  converged_x += sizeX-cutoff
  else
  update!(iterator.constr!, iterator.XBlocks.block, iterator.XBlocks.B_block)
- rand!(X)
- rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
+ rand!(rng, X)
+ rnext = lobpcg!(iterator, rng=rng, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
  append!(r, rnext, converged_x)
  converged_x += sizeX
  end

src/simple.jl

@@ -56,13 +56,14 @@ function powm_iterable!(A, x; tol = eps(real(eltype(A))) * size(A, 2) ^ 3, maxit
 end
 
 """
- powm(B; kwargs...) -> λ, x, [history]
+ powm(B; rng, kwargs...) -> λ, x, [history]
 
 See [`powm!`](@ref). Calls `powm!(B, x0; kwargs...)` with
-`x0` initialized as a random, complex unit vector.
+`x0` initialized as a random, complex unit vector using
+`rng::AbstractRNG`.
 """
-function powm(B; kwargs...)
- x0 = rand(Complex{real(eltype(B))}, size(B, 1))
+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...)
 end
@@ -149,13 +150,13 @@ function powm!(B, x;
 end
 
 """
- invpowm(B; shift = σ, kwargs...) -> λ, x, [history]
+ invpowm(B; shift = σ, rng, kwargs...) -> λ, x, [history]
 
 Find the approximate eigenpair `(λ, x)` of ``A`` near `shift`, where `B`
 is a linear map that has the effect `B * v = inv(A - σI) * v`.
 
 The method calls `powm!(B, x0; inverse = true, shift = σ)` with
-`x0` a random, complex unit vector. See [`powm!`](@ref)
+`x0` a random, complex unit vector using `rng::AbstractRNG`. See [`powm!`](@ref)
 
 # Examples
 
@@ -168,8 +169,8 @@ Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), 50, ismutating = true)
 λ, x = invpowm(Fmap, shift = σ, tol = 1e-4, maxiter = 200)
 ```
 """
-function invpowm(B; kwargs...)
- x0 = rand(Complex{real(eltype(B))}, size(B, 1))
+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...)
 end

src/svdl.jl

@@ -101,6 +101,7 @@ If `log` is set to `true` is given, method will output a tuple `X, L, ch`. Where
 
 - `nsv::Int = 6`: number of singular values requested;
 - `v0 = random unit vector`: starting guess vector in the domain of `A`.
+- `rng::AbstractRNG`: generator for pseudorandom initialization of `v0`
  The length of `q` should be the number of columns in `A`;
 - `k::Int = 2nsv`: maximum number of Lanczos vectors to compute before restarting;
 - `j::Int = nsv`: number of vectors to keep at the end of the restart.
@@ -175,7 +176,7 @@ end
 #########################
 
 function svdl_method!(log::ConvergenceHistory, A, l::Int=min(6, size(A,1)); k::Int=2l,
- j::Int=l, v0::AbstractVector = Vector{eltype(A)}(randn(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)