From cef2f42da57b58bfdf697ffe2708caaef6843422 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Sat, 30 Oct 2021 08:53:23 -0400 Subject: [PATCH 1/6] use deterministic RNG initial state --- Project.toml | 2 +- docs/src/eigenproblems/lobpcg.md | 4 ++++ docs/src/linear_systems/bicgstabl.md | 11 ++++++++--- docs/src/linear_systems/idrs.md | 7 ++++++- docs/src/svd/svdl.md | 6 ++++++ src/bicgstabl.jl | 6 ++++-- src/idrs.jl | 7 ++++--- src/lobpcg.jl | 19 ++++++++++++------- src/simple.jl | 17 +++++++++-------- src/svdl.jl | 3 ++- 10 files changed, 56 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index 17e24976..ee82ccf4 100644 --- a/Project.toml +++ b/Project.toml @@ -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.4" diff --git a/docs/src/eigenproblems/lobpcg.md b/docs/src/eigenproblems/lobpcg.md index 609ca0fd..b305b3cf 100644 --- a/docs/src/eigenproblems/lobpcg.md +++ b/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]. diff --git a/docs/src/linear_systems/bicgstabl.md b/docs/src/linear_systems/bicgstabl.md index 518a5317..d49ec5ab 100644 --- a/docs/src/linear_systems/bicgstabl.md +++ b/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. \ No newline at end of file diff --git a/docs/src/linear_systems/idrs.md b/docs/src/linear_systems/idrs.md index f530563d..106569e6 100644 --- a/docs/src/linear_systems/idrs.md +++ b/docs/src/linear_systems/idrs.md @@ -14,4 +14,9 @@ idrs! The current implementation is based on the [MATLAB version](http://ta.twi.tudelft.nl/nw/users/gijzen/idrs.m) by Van Gijzen and Sonneveld. For background see [^Sonneveld2008], [^VanGijzen2011] and [the IDR(s) webpage](http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html). [^Sonneveld2008]: IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems. P. Sonneveld and M. B. van Gijzen SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035--1062, 2008 -[^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 \ No newline at end of file +[^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 + +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. \ No newline at end of file diff --git a/docs/src/svd/svdl.md b/docs/src/svd/svdl.md index bb3915eb..f2faa87a 100644 --- a/docs/src/svd/svdl.md +++ b/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. diff --git a/src/bicgstabl.jl b/src/bicgstabl.jl index a13e330a..01df0104 100644 --- a/src/bicgstabl.jl +++ b/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(0), 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 diff --git a/src/idrs.jl b/src/idrs.jl index 1c827632..1efc8724 100644 --- a/src/idrs.jl +++ b/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(0)) 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) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 57eefec2..0e1afea5 100644 --- a/src/lobpcg.jl +++ b/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(0), 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(0), 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,13 +866,14 @@ end """ function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=false, + rng::AbstractRNG=MersenneTwister(0), 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) @@ -877,7 +881,7 @@ function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=fal 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(0), C=nothing, tol::Real=default_tolerance(eltype(X0))) n = size(X0, 1) sizeX = size(X0, 2) @@ -946,13 +951,13 @@ function lobpcg(A, B, largest::Bool, X0, nev::Int; 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]) + rand!(rng, X[:, cutoff+1:sizeX]) rnext = lobpcg!(iterator, 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) + rand!(rng, X) rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true) append!(r, rnext, converged_x) converged_x += sizeX diff --git a/src/simple.jl b/src/simple.jl index 326025c7..a5908289 100644 --- a/src/simple.jl +++ b/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(0), 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(0), kwargs...) + x0 = rand(rng, Complex{real(eltype(B))}, size(B, 1)) rmul!(x0, one(eltype(B)) / norm(x0)) invpowm!(B, x0; kwargs...) end diff --git a/src/svdl.jl b/src/svdl.jl index 908f4506..90b72006 100644 --- a/src/svdl.jl +++ b/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(0), 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) From c6f42de29554bc27d27efa7bec5ec3cd9b3811e1 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Sat, 30 Oct 2021 13:23:58 -0400 Subject: [PATCH 2/6] fix test failure --- src/svdl.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/svdl.jl b/src/svdl.jl index 90b72006..da988b00 100644 --- a/src/svdl.jl +++ b/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(0), v0::AbstractVector = Vector{eltype(A)}(randn(rng, size(A, 2))) |> x->rmul!(x, inv(norm(x))), + j::Int=l, rng::AbstractRNG=MersenneTwister(l), 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) From d89d1c804696aa2326c22de036d34291181b4d20 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Sat, 30 Oct 2021 16:36:48 -0400 Subject: [PATCH 3/6] use a better seed than 0 --- src/IterativeSolvers.jl | 4 ++++ src/bicgstabl.jl | 2 +- src/idrs.jl | 2 +- src/lobpcg.jl | 8 ++++---- src/simple.jl | 4 ++-- src/svdl.jl | 2 +- 6 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/IterativeSolvers.jl b/src/IterativeSolvers.jl index c3f11860..57a0f765 100644 --- a/src/IterativeSolvers.jl +++ b/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") diff --git a/src/bicgstabl.jl b/src/bicgstabl.jl index 01df0104..56f3b6b5 100644 --- a/src/bicgstabl.jl +++ b/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) diff --git a/src/idrs.jl b/src/idrs.jl index 1efc8724..6204ae63 100644 --- a/src/idrs.jl +++ b/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 diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 0e1afea5..525947e9 100644 --- a/src/lobpcg.jl +++ b/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) diff --git a/src/simple.jl b/src/simple.jl index a5908289..6ddc8bee 100644 --- a/src/simple.jl +++ b/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...) diff --git a/src/svdl.jl b/src/svdl.jl index da988b00..60fac7db 100644 --- a/src/svdl.jl +++ b/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) From d0013834f501b108ad8cae66f532f1a0f190e4f2 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Sat, 30 Oct 2021 16:50:49 -0400 Subject: [PATCH 4/6] pass rng recursively in lobpcg! --- src/lobpcg.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 525947e9..511c825e 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -943,7 +943,7 @@ 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 @@ -952,13 +952,13 @@ function lobpcg(A, B, largest::Bool, X0, nev::Int; update!(iterator.constr!, iterator.XBlocks.block[:, 1:cutoff], iterator.XBlocks.B_block[:, 1:cutoff]) X[:, 1:sizeX-cutoff] .= X[:, cutoff+1:sizeX] rand!(rng, X[:, cutoff+1:sizeX]) - rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true) + 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!(rng, X) - rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true) + rnext = lobpcg!(iterator, rng=rng, log=log, tol=tol, maxiter=maxiter, not_zeros=true) append!(r, rnext, converged_x) converged_x += sizeX end From 5af2255683ef75304261115f6fdabc75010808f5 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 14 Dec 2021 17:23:25 -0500 Subject: [PATCH 5/6] Update Project.toml Co-authored-by: Moritz Schauer --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ee82ccf4..f2e10848 100644 --- a/Project.toml +++ b/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" From 5f679713c39b453ae8efbf54b58f089fef7d7cf3 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 29 Jun 2022 08:51:40 +0200 Subject: [PATCH 6/6] Require Julia 1.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f2e10848..083af3c2 100644 --- a/Project.toml +++ b/Project.toml @@ -11,4 +11,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] RecipesBase = "0.6, 0.7, 0.8, 1.0" -julia = "1.4" +julia = "1.6"