From e9cabd35d4b9a24391c4d72ee515e4d136f43191 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 26 Oct 2015 13:36:10 -0400 Subject: [PATCH 1/2] Fix promotion for eigs and svds --- base/linalg/arnoldi.jl | 43 ++++++++++++++++++++++++++++++++++-------- test/linalg/arnoldi.jl | 8 ++++++++ 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/base/linalg/arnoldi.jl b/base/linalg/arnoldi.jl index a6c843605eed1..4d3c74554889d 100644 --- a/base/linalg/arnoldi.jl +++ b/base/linalg/arnoldi.jl @@ -48,8 +48,20 @@ The following keyword arguments are supported: ``` """ eigs(A; args...) = eigs(A, I; args...) - - +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) @@ -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) @@ -229,12 +242,19 @@ 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 ] @@ -242,7 +262,14 @@ end 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 @@ -250,7 +277,7 @@ function svds{S}(X::S; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, m 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]) diff --git a/test/linalg/arnoldi.jl b/test/linalg/arnoldi.jl index 78b14e84b2ba8..89382dbb52e4c 100644 --- a/test/linalg/arnoldi.jl +++ b/test/linalg/arnoldi.jl @@ -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))) From fa006bcd3dad533ab1b789fed6e7a042b172c622 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 26 Oct 2015 14:30:09 -0400 Subject: [PATCH 2/2] Improve error messages from ARPACK. Fixes #13709 --- base/linalg/arnoldi.jl | 32 ++++++++++++++++---------------- base/linalg/arpack.jl | 22 ++++++++++++++-------- base/linalg/exceptions.jl | 13 ++++++++++++- 3 files changed, 42 insertions(+), 25 deletions(-) diff --git a/base/linalg/arnoldi.jl b/base/linalg/arnoldi.jl index 4d3c74554889d..2eaf5bbc76c35 100644 --- a/base/linalg/arnoldi.jl +++ b/base/linalg/arnoldi.jl @@ -47,20 +47,20 @@ The following keyword arguments are supported: =============== ================================== ================================== ``` """ -eigs(A; args...) = eigs(A, I; args...) -eigs{T<:BlasFloat}(A::AbstractMatrix{T}, ::UniformScaling; 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}; 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...) +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; args...) + eigs(convert(AbstractMatrix{Tnew}, A), I; kwargs...) end -function eigs(A::AbstractMatrix, B::AbstractMatrix; args...) +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); args...) + eigs(convert(AbstractMatrix{Tnew}, A), convert(AbstractMatrix{Tnew}, B); kwargs...) end doc""" ```rst @@ -106,7 +106,7 @@ The following keyword arguments are supported: =============== ================================== ================================== ``` """ -eigs(A, B; args...) = _eigs(A, B; args...) +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,)), @@ -147,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")) @@ -262,13 +262,13 @@ end size(s::SVDOperator) = s.m + s.n, s.m + s.n issym(s::SVDOperator) = true -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...) +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); args...) + svds(convert(AbstractMatrix{Tnew}, A); kwargs...) end -svds(A; args...) = _svds(A; args...) +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")) diff --git a/base/linalg/arpack.jl b/base/linalg/arpack.jl index 9b31ab38ecb9a..ce564e2553b7e 100644 --- a/base/linalg/arpack.jl +++ b/base/linalg/arpack.jl @@ -59,7 +59,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 == "I" # corresponds to dsdrv2, dndrv2 or zndrv2 if ido[1] == -1 || ido[1] == 1 @@ -67,7 +67,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 == 2 # corresponds to dsdrv3, dndrv3 or zndrv3 if ido[1] == -1 || ido[1] == 1 @@ -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 @@ -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")) @@ -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) @@ -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) @@ -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 @@ -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 diff --git a/base/linalg/exceptions.jl b/base/linalg/exceptions.jl index bd7686193c11e..c8975e45d53fc 100644 --- a/base/linalg/exceptions.jl +++ b/base/linalg/exceptions.jl @@ -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