Skip to content

Commit

Permalink
Fix promotion for eigs and svds
Browse files Browse the repository at this point in the history
  • Loading branch information
andreasnoack committed Oct 26, 2015
1 parent 77b2527 commit e9cabd3
Show file tree
Hide file tree
Showing 2 changed files with 43 additions and 8 deletions.
43 changes: 35 additions & 8 deletions base/linalg/arnoldi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,20 @@ The following keyword arguments are supported:
```
"""
eigs(A; args...) = eigs(A, I; args...)

This comment has been minimized.

Copy link
@jiahao

jiahao Oct 26, 2015

Member

By convention, keyword arguments are splatted as ; kwargs... everywhere else in Base except here. Could you rename args to kwargs in these methods?



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

eigs{T<:BlasFloat}(A::AbstractMatrix{T}, B::AbstractMatrix{T}; args...) = _eigs(A, B; args...)
eigs(A::AbstractMatrix{BigFloat}, B::AbstractMatrix...; args...) = throw(MethodError(eigs, Any[A,B,args...]))
eigs(A::AbstractMatrix{BigFloat}, B::UniformScaling; args...) = throw(MethodError(eigs, Any[A,B,args...]))
function eigs{T}(A::AbstractMatrix{T}, ::UniformScaling; args...)
Tnew = typeof(zero(T)/sqrt(one(T)))
eigs(convert(AbstractMatrix{Tnew}, A), I; args...)
end
function eigs(A::AbstractMatrix, B::AbstractMatrix; args...)
T = promote_type(eltype(A), eltype(B))
Tnew = typeof(zero(T)/sqrt(one(T)))
eigs(convert(AbstractMatrix{Tnew}, A), convert(AbstractMatrix{Tnew}, B); args...)
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; args...) = _eigs(A, B; args...)
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 @@ -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}; args...) = _svds(A; args...)
svds(A::AbstractMatrix{BigFloat}; args...) = throw(MethodError(svds, Any[A, args...]))
function svds{T}(A::AbstractMatrix{T}; args...)
Tnew = typeof(zero(T)/sqrt(one(T)))
svds(convert(AbstractMatrix{Tnew}, A); args...)
end
svds(A; args...) = _svds(A; args...)
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
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 e9cabd3

Please sign in to comment.