Skip to content

Commit

Permalink
Tolerance increase, test matrix change, and system test tolerance
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohamed82008
mohamed82008 committed Jul 12, 2018
1 parent 750b7ff commit 42a435f
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 66 deletions.
29 changes: 14 additions & 15 deletions src/lobpcg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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")
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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.")
Expand Down
97 changes: 46 additions & 51 deletions test/lobpcg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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

Expand Down

0 comments on commit 42a435f

Please sign in to comment.