Skip to content

Commit

Permalink
Merge pull request #13778 from JuliaLang/anj/arpack
Browse files Browse the repository at this point in the history
Some promotion and error reporting fixes for eigs
  • Loading branch information
andreasnoack committed Oct 28, 2015
2 parents f2ac794 + fa006bc commit 361cf34
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 19 deletions.
47 changes: 37 additions & 10 deletions base/linalg/arnoldi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,21 @@ The following keyword arguments are supported:
=============== ================================== ==================================
```
"""
eigs(A; args...) = eigs(A, I; args...)


eigs(A; kwargs...) = eigs(A, I; kwargs...)
eigs{T<:BlasFloat}(A::AbstractMatrix{T}, ::UniformScaling; kwargs...) = _eigs(A, I; kwargs...)

eigs{T<:BlasFloat}(A::AbstractMatrix{T}, B::AbstractMatrix{T}; kwargs...) = _eigs(A, B; kwargs...)
eigs(A::AbstractMatrix{BigFloat}, B::AbstractMatrix...; kwargs...) = throw(MethodError(eigs, Any[A,B,kwargs...]))
eigs(A::AbstractMatrix{BigFloat}, B::UniformScaling; kwargs...) = throw(MethodError(eigs, Any[A,B,kwargs...]))
function eigs{T}(A::AbstractMatrix{T}, ::UniformScaling; kwargs...)
Tnew = typeof(zero(T)/sqrt(one(T)))
eigs(convert(AbstractMatrix{Tnew}, A), I; kwargs...)
end
function eigs(A::AbstractMatrix, B::AbstractMatrix; kwargs...)
T = promote_type(eltype(A), eltype(B))
Tnew = typeof(zero(T)/sqrt(one(T)))
eigs(convert(AbstractMatrix{Tnew}, A), convert(AbstractMatrix{Tnew}, B); kwargs...)
end
doc"""
```rst
.. eigs(A, B; nev=6, ncv=max(20,2*nev+1), which="LM", tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)
Expand Down Expand Up @@ -94,7 +106,8 @@ The following keyword arguments are supported:
=============== ================================== ==================================
```
"""
function eigs(A, B;
eigs(A, B; kwargs...) = _eigs(A, B; kwargs...)
function _eigs(A, B;
nev::Integer=6, ncv::Integer=max(20,2*nev+1), which=:LM,
tol=0.0, maxiter::Integer=300, sigma=nothing, v0::Vector=zeros(eltype(A),(0,)),
ritzvec::Bool=true)
Expand Down Expand Up @@ -134,7 +147,7 @@ function eigs(A, B;
end
if (which != :LM && which != :SM && which != :LR && which != :SR &&
which != :LI && which != :SI && which != :BE)
throw(ArgumentError("which must be :LM, :SM, :LR, :SR, :LI, :SI, or :BE, got $(repr(which))"))
throw(ArgumentError("which must be :LM, :SM, :LR, :SR, :LI, :SI, or :BE, got $(repr(which))"))
end
if which == :BE && !sym
throw(ArgumentError("which=:BE only possible for real symmetric problem"))
Expand Down Expand Up @@ -229,28 +242,42 @@ end


## svds

type SVDOperator{T,S} <: AbstractArray{T, 2}
### Restrict operator to BlasFloat because ARPACK only supports that. Loosen restriction
### when we switch to our own implementation
type SVDOperator{T<:BlasFloat,S} <: AbstractArray{T, 2}
X::S
m::Int
n::Int
SVDOperator(X::S) = new(X, size(X,1), size(X,2))
SVDOperator(X::AbstractMatrix) = new(X, size(X, 1), size(X, 2))
end

function SVDOperator{T}(A::AbstractMatrix{T})
Tnew = typeof(zero(T)/sqrt(one(T)))
Anew = convert(AbstractMatrix{Tnew}, A)
SVDOperator{Tnew,typeof(Anew)}(Anew)
end

## v = [ left_singular_vector; right_singular_vector ]
*{T,S}(s::SVDOperator{T,S}, v::Vector{T}) = [s.X * v[s.m+1:end]; s.X' * v[1:s.m]]
size(s::SVDOperator) = s.m + s.n, s.m + s.n
issym(s::SVDOperator) = true

function svds{S}(X::S; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxiter::Int = 1000)
svds{T<:BlasFloat}(A::AbstractMatrix{T}; kwargs...) = _svds(A; kwargs...)
svds(A::AbstractMatrix{BigFloat}; kwargs...) = throw(MethodError(svds, Any[A, kwargs...]))
function svds{T}(A::AbstractMatrix{T}; kwargs...)
Tnew = typeof(zero(T)/sqrt(one(T)))
svds(convert(AbstractMatrix{Tnew}, A); kwargs...)
end
svds(A; kwargs...) = _svds(A; kwargs...)
function _svds(X; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxiter::Int = 1000)
if nsv < 1
throw(ArgumentError("number of singular values (nsv) must be ≥ 1, got $nsv"))
end
if nsv > minimum(size(X))
throw(ArgumentError("number of singular values (nsv) must be ≤ $(minimum(size(X))), got $nsv"))
end
otype = eltype(X)
ex = eigs(SVDOperator{otype,S}(X), I; ritzvec = ritzvec, nev = 2*nsv, tol = tol, maxiter = maxiter)
ex = eigs(SVDOperator(X), I; ritzvec = ritzvec, nev = 2*nsv, tol = tol, maxiter = maxiter)
ind = [1:2:nsv*2;]
sval = abs(ex[1][ind])

Expand Down
22 changes: 14 additions & 8 deletions base/linalg/arpack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -59,15 +59,15 @@ function aupd_wrapper(T, matvecA::Function, matvecB::Function, solveSI::Function
elseif ido[1] == 99
break
else
error("Internal ARPACK error")
throw(ARPACKException("unexpected behavior"))
end
elseif mode == 3 && bmat == "I" # corresponds to dsdrv2, dndrv2 or zndrv2
if ido[1] == -1 || ido[1] == 1
workd[store_idx] = solveSI(x)
elseif ido[1] == 99
break
else
error("Internal ARPACK error")
throw(ARPACKException("unexpected behavior"))
end
elseif mode == 2 # corresponds to dsdrv3, dndrv3 or zndrv3
if ido[1] == -1 || ido[1] == 1
Expand All @@ -81,7 +81,7 @@ function aupd_wrapper(T, matvecA::Function, matvecB::Function, solveSI::Function
elseif ido[1] == 99
break
else
error("Internal ARPACK error")
throw(ARPACKException("unexpected behavior"))
end
elseif mode == 3 && bmat == "G" # corresponds to dsdrv4, dndrv4 or zndrv4
if ido[1] == -1
Expand All @@ -93,7 +93,7 @@ function aupd_wrapper(T, matvecA::Function, matvecB::Function, solveSI::Function
elseif ido[1] == 99
break
else
error("Internal ARPACK error")
throw(ARPACKException("unexpected behavior"))
end
else
throw(ArgumentError("ARPACK mode ($mode) not yet supported"))
Expand Down Expand Up @@ -133,7 +133,9 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ByteString,
neupd(ritzvec, howmny, select, d, v, ldv, sigmar, workev,
bmat, n, which, nev, TOL, resid, ncv, v, ldv,
iparam, ipntr, workd, workl, lworkl, rwork, info)
if info[1] != 0; throw(ARPACKException(info[1])); end
if info[1] != 0
throw(ARPACKException(info[1]))
end

p = sortperm(dmap(d[1:nev]), rev=true)
return ritzvec ? (d[p], v[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d[p],iparam[5],iparam[3],iparam[9],resid)
Expand All @@ -145,7 +147,9 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ByteString,
seupd(ritzvec, howmny, select, d, v, ldv, sigmar,
bmat, n, which, nev, TOL, resid, ncv, v, ldv,
iparam, ipntr, workd, workl, lworkl, info)
if info[1] != 0; throw(ARPACKException(info[1])); end
if info[1] != 0
throw(ARPACKException(info[1]))
end

p = sortperm(dmap(d), rev=true)
return ritzvec ? (d[p], v[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d,iparam[5],iparam[3],iparam[9],resid)
Expand All @@ -162,7 +166,9 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ByteString,
neupd(ritzvec, howmny, select, dr, di, v, ldv, sigmar, sigmai,
workev, bmat, n, which, nev, TOL, resid, ncv, v, ldv,
iparam, ipntr, workd, workl, lworkl, info)
if info[1] != 0; throw(ARPACKException(info[1])); end
if info[1] != 0
throw(ARPACKException(info[1]))
end
evec = complex(Array(T, n, nev+1), Array(T, n, nev+1))

j = 1
Expand All @@ -181,7 +187,7 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ByteString,
evec[:,j] = v[:,j]
j += 1
else
error("ARPACK unexpected behavior")
throw(ARPACKException("unexpected behavior"))
end
end

Expand Down
13 changes: 12 additions & 1 deletion base/linalg/exceptions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,18 @@ type LAPACKException <: Exception
end

type ARPACKException <: Exception
info::BlasInt
info::ByteString
end

function ARPACKException(i::Integer)
if i == -8
return ARPACKException("error return from calculation of a real Schur form.")
elseif i == -9
return ARPACKException("error return from calculation of eigenvectors.")
elseif i == -14
return ARPACKException("did not find any eigenvalues to sufficient accuracy. Try with a different starting vector or more Lanczos vectors by increasing the value of ncv.")
end
return ARPACKException("unspecified ARPACK error: $i")
end

type SingularException <: Exception
Expand Down
8 changes: 8 additions & 0 deletions test/linalg/arnoldi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -204,3 +204,11 @@ let # complex svds test
@test_throws ArgumentError svds(A,nsv=0)
@test_throws ArgumentError svds(A,nsv=20)
end

# test promotion
eigs(rand(1:10, 10, 10))
eigs(rand(1:10, 10, 10), rand(1:10, 10, 10) |> t -> t't)
svds(rand(1:10, 10, 8))
@test_throws MethodError eigs(big(rand(1:10, 10, 10)))
@test_throws MethodError eigs(big(rand(1:10, 10, 10)), rand(1:10, 10, 10))
@test_throws MethodError svds(big(rand(1:10, 10, 8)))

0 comments on commit 361cf34

Please sign in to comment.