`DomainError` in `svdsolve` when running multithreaded code

[alanderos91]

Sometimes svdsolve throws a DomainError on v0.5.4 and now v0.6.0:

using LinearAlgebra, KrylovKit, StableRNGs

m, n = 4100, 4200
X = randn(StableRNG(1234), m, n)

# should work for r <= 4096 == BLOCKSIZE
r = 10
svdsolve(X, m, r, :LR, Float64, krylovdim=r);

# this should throw the DomainError; takes a long time
r = 4097 # BLOCKSIZE + 1
svdsolve(X, m, r, :LR, Float64, krylovdim=r);

The stacktrace (see below) hints at a prevpow call that may be the source of the problem; see for example https://github.com/Jutho/KrylovKit.jl/blob/v0.6.0/src/orthonormal.jl#L167. My guess is that we somehow end up with prevpow(2, 0) which throws the error.

I did not study the multithreaded code closely but would adding something like

prevpow(2, max(1, div(BLOCKSIZE, n)))

be sufficient? It certainly avoids the error but not sure about impact on performance. Admittedly, this is a highly contrived edge case. I can't recall how I came across this error in the first place a few months ago. My example was constructed after peeking at the multithreaded code.

Other information

Stacktrace:

ERROR: DomainError with 0:
`x` must be ≥ 1.
Stacktrace:
  [1] prevpow(a::Int64, x::Int64)
    @ Base ./intfuncs.jl:479
  [2] unproject_linear_multithreaded!(y::Vector{Float64}, b::KrylovKit.OrthonormalBasis{Vector{Float64}}, x::SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, α::Int64, β::Int64, r::Base.OneTo{Int64})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/orthonormal.jl:167
  [3] unproject!(y::Vector{Float64}, b::KrylovKit.OrthonormalBasis{Vector{Float64}}, x::SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, α::Int64, β::Int64, r::Base.OneTo{Int64})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/orthonormal.jl:135
  [4] unproject!
    @ ~/.julia/packages/KrylovKit/diNbc/src/orthonormal.jl:134 [inlined]
  [5] mul!
    @ ~/.julia/packages/KrylovKit/diNbc/src/orthonormal.jl:61 [inlined]
  [6] *
    @ ~/.julia/packages/KrylovKit/diNbc/src/orthonormal.jl:59 [inlined]
  [7] #69
    @ ./none:0 [inlined]
  [8] iterate
    @ ./generator.jl:47 [inlined]
  [9] collect(itr::Base.Generator{KrylovKit.ColumnIterator{SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, Base.OneTo{Int64}}, KrylovKit.var"#69#73"{KrylovKit.OrthonormalBasis{Vector{Float64}}}})
    @ Base ./array.jl:724
 [10] svdsolve(A::Matrix{Float64}, x₀::Vector{Float64}, howmany::Int64, which::Symbol, alg::GKL{ModifiedGramSchmidt2, Float64})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/svdsolve.jl:274
 [11] #svdsolve#68
    @ ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/svdsolve.jl:127 [inlined]
 [12] svdsolve(f::Matrix{Float64}, n::Int64, howmany::Int64, which::Symbol, T::Type; kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:krylovdim,), Tuple{Int64}}})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/svdsolve.jl:119
 [13] top-level scope
    @ REPL[16]:1

Julia command:

julia -t 10

Version info:

julia> versioninfo()
Julia Version 1.7.3
Commit 742b9abb4d (2022-05-06 12:58 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i9-10900KF CPU @ 3.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

Environment:

(@v1.7) pkg> activate --temp
  Activating new project at `/tmp/jl_NkV4GT`

(jl_NkV4GT) pkg> add LinearAlgebra KrylovKit StableRNGs
   Resolving package versions...
    Updating `/tmp/jl_NkV4GT/Project.toml`
  [0b1a1467] + KrylovKit v0.6.0
  [860ef19b] + StableRNGs v1.0.0
  [37e2e46d] + LinearAlgebra
    Updating `/tmp/jl_NkV4GT/Manifest.toml`
  [79e6a3ab] + Adapt v3.5.0
  [d360d2e6] + ChainRulesCore v1.15.7
  [34da2185] + Compat v4.6.0
  [46192b85] + GPUArraysCore v0.1.3
  [0b1a1467] + KrylovKit v0.6.0
  [860ef19b] + StableRNGs v1.0.0
  [56f22d72] + Artifacts
  [2a0f44e3] + Base64
  [ade2ca70] + Dates
  [b77e0a4c] + InteractiveUtils
  [8f399da3] + Libdl
  [37e2e46d] + LinearAlgebra
  [56ddb016] + Logging
  [d6f4376e] + Markdown
  [de0858da] + Printf
  [9a3f8284] + Random
  [ea8e919c] + SHA
  [9e88b42a] + Serialization
  [2f01184e] + SparseArrays
  [8dfed614] + Test
  [cf7118a7] + UUIDs
  [4ec0a83e] + Unicode
  [e66e0078] + CompilerSupportLibraries_jll
  [4536629a] + OpenBLAS_jll
  [8e850b90] + libblastrampoline_jll