diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 848c73b3..b5009e37 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -743,6 +743,8 @@ function (iterator::LOBPCGIterator{Generalized})(residualTolerance, log) where { end end +default_tolerance(::Type{T}) where {T<:Number} = eps(real(T))^(real(T)(3)/10) + """ The Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG) @@ -771,7 +773,7 @@ Finds the `nev` extremal eigenvalues and their corresponding eigenvectors satisf - `maxiter`: maximum number of iterations; default is 200; -- `tol::Number`: tolerance to which residual vector norms must be under. +- `tol::Real`: tolerance to which residual vector norms must be under. # Output @@ -808,7 +810,7 @@ end - `maxiter`: maximum number of iterations; default is 200; -- `tol::Number`: tolerance to which residual vector norms must be under. +- `tol::Real`: tolerance to which residual vector norms must be under. # Output @@ -818,10 +820,9 @@ function lobpcg(A, largest::Bool, X0; kwargs...) lobpcg(A, nothing, largest, X0; kwargs...) end function lobpcg(A, B, largest, X0; - not_zeros=false, log=false, P=nothing, - C=nothing, tol=nothing, maxiter=200) + not_zeros=false, log=false, P=nothing, maxiter=200, + C=nothing, tol::Real=default_tolerance(eltype(X0))) X = copy(X0) - T = eltype(X) n = size(X, 1) sizeX = size(X, 2) sizeX > n && throw("X column dimension exceeds the row dimension") @@ -848,15 +849,15 @@ end - `maxiter`: maximum number of iterations; default is 200; -- `tol::Number`: tolerance to which residual vector norms must be under. +- `tol::Real`: tolerance to which residual vector norms must be under. # Output - `results`: a `LOBPCGResults` struct. `r.λ` and `r.X` store the eigenvalues and eigenvectors. """ -function lobpcg!(iterator::LOBPCGIterator; log=false, tol=nothing, maxiter=200, not_zeros=false) - T = eltype(iterator.XBlocks.block) +function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=false, + tol::Real=default_tolerance(eltype(iterator.XBlocks.block))) X = iterator.XBlocks.block iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block) if !not_zeros @@ -869,10 +870,9 @@ function lobpcg!(iterator::LOBPCGIterator; log=false, tol=nothing, maxiter=200, end n = size(X, 1) sizeX = size(X, 2) - residualTolerance = (tol isa Void) ? (eps(real(T)))^(real(T)(4)/10) : real(tol) iterator.iteration[] = 1 while iterator.iteration[] <= maxiter - state = iterator(residualTolerance, log) + state = iterator(tol, log) if log push!(iterator.trace, state) end @@ -881,7 +881,7 @@ function lobpcg!(iterator::LOBPCGIterator; log=false, tol=nothing, maxiter=200, end @inbounds iterator.λ .= view(iterator.ritz_values, 1:sizeX) - results = LOBPCGResults(iterator.λ, X, residualTolerance, iterator.residuals, iterator.iteration[], maxiter, all((x)->(norm(x)<=residualTolerance), view(iterator.residuals, 1:sizeX)), iterator.trace) + results = LOBPCGResults(iterator.λ, X, tol, iterator.residuals, iterator.iteration[], maxiter, all((x)->(norm(x)<=tol), view(iterator.residuals, 1:sizeX)), iterator.trace) return results end @@ -910,7 +910,7 @@ lobpcg(A, [B,] largest, X0, nev; kwargs...) -> results - `maxiter`: maximum number of iterations; default is 200; -- `tol::Number`: tolerance to which residual vector norms must be under. +- `tol::Real`: tolerance to which residual vector norms must be under. # Output @@ -920,9 +920,8 @@ function lobpcg(A, largest::Bool, X0, nev::Int; kwargs...) lobpcg(A, nothing, largest, X0, nev; kwargs...) end function lobpcg(A, B, largest::Bool, X0, nev::Int; - not_zeros=false, log=false, P=nothing, - C=nothing, tol=nothing, maxiter=200) - T = eltype(X0) + not_zeros=false, log=false, P=nothing, maxiter=200, + C=nothing, tol::Real=default_tolerance(eltype(X0))) n = size(X0, 1) sizeX = size(X0, 2) nev > n && throw("Number of eigenvectors desired exceeds the row dimension.") diff --git a/test/lobpcg.jl b/test/lobpcg.jl index fab7cfc0..c5e3bd29 100644 --- a/test/lobpcg.jl +++ b/test/lobpcg.jl @@ -33,10 +33,9 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I b = rand(T, n, 1) - tol = √eps(real(T)) - + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, largest, b; tol=tol, maxiter=Inf, log=false) λ, X = r.λ, r.X @test norm(A*X - X*λ) ≤ tol @@ -51,12 +50,11 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I + B = B' + B + 20I b = rand(T, n, 1) - tol = √eps(real(T)) - + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, b; tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @test max_err(A*X - B*X*λ) ≤ tol @@ -72,8 +70,7 @@ end A = laplace_matrix(Float64, 20, 2) rhs = randn(size(A, 2), 1) scale!(rhs, inv(norm(rhs))) - tol = 1e-5 - + tol = IterativeSolvers.default_tolerance(Float64) @testset "Matrix" begin @testset "largest = $largest" for largest in (true, false) r = lobpcg(A, largest, rhs; tol=tol, maxiter=Inf) @@ -87,10 +84,9 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I b = zeros(T, n, 1) - tol = √eps(real(T)) - + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, largest, b; tol=tol, maxiter=Inf, log=false) λ, X = r.λ, r.X @test norm(A*X - X*λ) ≤ tol @@ -101,11 +97,11 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I + B = B' + B + 20I b = zeros(T, n, 1) - tol = √eps(real(T)) + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, b; tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @@ -119,8 +115,8 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I - tol = √eps(real(T)) + A = A' + A + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false) λ, X = r.λ, r.X @@ -132,10 +128,10 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I - tol = √eps(real(T)) + B = B' + B + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @@ -149,8 +145,8 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I - tol = √eps(real(T)) + A = A' + A + 20I + tol = IterativeSolvers.default_tolerance(T) b = rand(T, n, 1) itr = LOBPCGIterator(A, largest, b) @@ -164,11 +160,11 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I + B = B' + B + 20I b = rand(T, n, 1) - tol = √eps(real(T)) + tol = IterativeSolvers.default_tolerance(T) itr = LOBPCGIterator(A, B, largest, b) r = lobpcg!(itr; tol=tol, maxiter=Inf, log=true) @@ -183,8 +179,8 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I - tol = √eps(real(T)) + A = A' + A + 20I + tol = IterativeSolvers.default_tolerance(T) P = JacobiPrec(diag(A)) r = lobpcg(A, largest, 1; P=P, tol=tol, maxiter=Inf, log=false) λ, X = r.λ, r.X @@ -196,11 +192,11 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I P = JacobiPrec(diag(A)) B = rand(T, n, n) - B = B' * B + I - tol = √eps(real(T)) + B = B' + B + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, 1; P=P, tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @@ -214,8 +210,8 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I - tol = √eps(real(T)) + A = A' + A + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false) λ1, X1 = r.λ, r.X r = lobpcg(A, largest, 1; C=copy(r.X), tol=tol, maxiter=Inf, log=false) @@ -229,10 +225,10 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I - tol = eps(real(T))^0.4 + B = B' + B + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=false) λ1, X1 = r.λ, r.X r = lobpcg(A, B, largest, 1; C=copy(r.X), tol=tol, maxiter=Inf, log=false) @@ -251,9 +247,9 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I b = rand(T, n, 2) - tol = √eps(real(T)) + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, largest, b; tol=tol, maxiter=Inf, log=false) λ, X = r.λ, r.X @@ -269,11 +265,11 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I + B = B' + B + 20I b = rand(T, n, 2) - tol = eps(real(T))^(real(T)(4/10)) + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, b; tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @test max_err(A*X - B*X*diagm(λ)) ≤ tol @@ -292,8 +288,8 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I - tol = eps(real(T))^0.4 + A = A' + A + 20I + tol = IterativeSolvers.default_tolerance(T) X0 = rand(T, n, block_size) r = lobpcg(A, largest, X0, 3, tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @@ -306,11 +302,10 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + I - tol = eps(real(T))^0.4 - + B = B' + B + 20I + tol = IterativeSolvers.default_tolerance(T) X0 = rand(T, n, block_size) r = lobpcg(A, B, largest, X0, 3, tol=tol, maxiter=Inf, log=true) λ, X = r.λ, r.X @@ -324,8 +319,8 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + I - tol = √eps(real(T)) + A = A' + A + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false) λ1, X1 = r.λ, r.X @@ -342,10 +337,10 @@ end @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) @testset "largest = $largest" for largest in (true, false) A = rand(T, n, n) - A = A' * A + 2I + A = A' + A + 20I B = rand(T, n, n) - B = B' * B + 2I - tol = eps(real(T))^0.4 + B = B' + B + 20I + tol = IterativeSolvers.default_tolerance(T) r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=false) λ1, X1 = r.λ, r.X