Skip to content

Commit

Permalink
Merge pull request #15549 from JuliaLang/nl/lapack3.6
Browse files Browse the repository at this point in the history 
Backport LAPACK 3.6 support
  • Loading branch information
@nalimilan
nalimilan committed Mar 19, 2016
2 parents a7d8a3a + e3f3f72 commit 03c072d
Show file tree
Hide file tree
Showing 4 changed files with 198 additions and 6 deletions.
165 changes: 163 additions & 2 deletions base/linalg/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,12 @@ import Base.blasfunc
import ..LinAlg: BlasFloat, Char, BlasInt, LAPACKException,
DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare

const VERSION = Ref{VersionNumber}()

function __init__()
VERSION[] = VersionNumber(laver()...)
end

#Generic LAPACK error handlers
macro assertargsok() #Handle only negative info codes - use only if positive info code is useful!
:(info[1]<0 && throw(ArgumentError("invalid argument #$(-info[1]) to LAPACK call")))
Expand Down Expand Up @@ -52,6 +58,26 @@ end
subsetrows(X::AbstractVector, Y::AbstractArray, k) = Y[1:k]
subsetrows(X::AbstractMatrix, Y::AbstractArray, k) = Y[1:k, :]

function checkfinite(A::StridedMatrix)
for i = eachindex(A)
if !isfinite(A[i])
throw(ArgumentError("matrix contains NaNs"))
end
end
return true
end

# LAPACK version number
@eval function laver()
major = BlasInt[0]
minor = BlasInt[0]
patch = BlasInt[0]
ccall(($(blasfunc(:ilaver_)), liblapack), Void,
(Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
major, minor, patch)
return major[1], minor[1], patch[1]
end

# (GB) general banded matrices, LU decomposition and solver
for (gbtrf, gbtrs, elty) in
((:dgbtrf_,:dgbtrs_,:Float64),
Expand Down Expand Up @@ -146,6 +172,7 @@ for (gebal, gebak, elty, relty) in
# DOUBLE PRECISION A( LDA, * ), SCALE( * )
function gebal!(job::Char, A::StridedMatrix{$elty})
chkstride1(A)
checkfinite(A) # balancing routines don't support NaNs and Infs
n = chksquare(A)
info = Array(BlasInt, 1)
ihi = Array(BlasInt, 1)
Expand All @@ -169,6 +196,7 @@ for (gebal, gebak, elty, relty) in
ilo::BlasInt, ihi::BlasInt, scale::Vector{$relty},
V::StridedMatrix{$elty})
chkstride1(V)
checkfinite(V) # balancing routines don't support NaNs and Infs
chkside(side)
n = chksquare(V)
info = Array(BlasInt, 1)
Expand Down Expand Up @@ -1382,6 +1410,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in
# $ WI( * ), WORK( * ), WR( * )
function geev!(jobvl::Char, jobvr::Char, A::StridedMatrix{$elty})
chkstride1(A)
checkfinite(A) # balancing routines don't support NaNs and Infs
n = chksquare(A)
lvecs = jobvl == 'V'
rvecs = jobvr == 'V'
Expand Down Expand Up @@ -1654,11 +1683,140 @@ Finds the generalized singular value decomposition of `A` and `B`, `U'*A*Q = D1*
and `V'*B*Q = D2*R`. `D1` has `alpha` on its diagonal and `D2` has `beta` on its
diagonal. If `jobu = U`, the orthogonal/unitary matrix `U` is computed. If
`jobv = V` the orthogonal/unitary matrix `V` is computed. If `jobq = Q`,
the orthogonal/unitary matrix `Q` is computed. If `job{u,v,q} = N`, that
matrix is not computed.
the orthogonal/unitary matrix `Q` is computed. If `jobu`, `jobv` or `jobq` is
`N`, that matrix is not computed. This function is only available in LAPACK
versions prior to 3.6.0.
"""
ggsvd!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix, B::StridedMatrix)


for (f, elty) in ((:dggsvd3_, :Float64),
(:sggsvd3_, :Float32))
@eval begin
function ggsvd3!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
chkstride1(A, B)
m, n = size(A)
if size(B, 2) != n
throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs $n"))
end
p = size(B, 1)
k = Ref{BlasInt}()
l = Ref{BlasInt}()
lda = max(1, stride(A, 2))
ldb = max(1, stride(B, 2))
alpha = similar(A, $elty, n)
beta = similar(A, $elty, n)
ldu = max(1, m)
U = jobu == 'U' ? similar(A, $elty, ldu, m) : similar(A, $elty, 0)
ldv = max(1, p)
V = jobv == 'V' ? similar(A, $elty, ldv, p) : similar(A, $elty, 0)
ldq = max(1, n)
Q = jobq == 'Q' ? similar(A, $elty, ldq, n) : similar(A, $elty, 0)
work = Array($elty, 1)
lwork = BlasInt(-1)
iwork = Array(BlasInt, n)
info = Ref{BlasInt}()
for i = 1:2
ccall(($(blasfunc(f)), liblapack), Void,
(Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}),
&jobu, &jobv, &jobq, &m,
&n, &p, k, l,
A, &lda, B, &ldb,
alpha, beta, U, &ldu,
V, &ldv, Q, &ldq,
work, &lwork, iwork, info)
@lapackerror
if i == 1
lwork = BlasInt(work[1])
resize!(work, lwork)
end
end
if m - k[] - l[] >= 0
R = triu(A[1:k[] + l[],n - k[] - l[] + 1:n])
else
R = triu([A[1:m, n - k[] - l[] + 1:n]; B[m - k[] + 1:l[], n - k[] - l[] + 1:n]])
end
return U, V, Q, alpha, beta, k[], l[], R
end
end
end

for (f, elty, relty) in ((:zggsvd3_, :Complex128, :Float64),
(:cggsvd3_, :Complex64, :Float32))
@eval begin
function ggsvd3!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
chkstride1(A, B)
m, n = size(A)
if size(B, 2) != n
throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs $n"))
end
p = size(B, 1)
k = Array(BlasInt, 1)
l = Array(BlasInt, 1)
lda = max(1,stride(A, 2))
ldb = max(1,stride(B, 2))
alpha = similar(A, $relty, n)
beta = similar(A, $relty, n)
ldu = max(1, m)
U = jobu == 'U' ? similar(A, $elty, ldu, m) : similar(A, $elty, 0)
ldv = max(1, p)
V = jobv == 'V' ? similar(A, $elty, ldv, p) : similar(A, $elty, 0)
ldq = max(1, n)
Q = jobq == 'Q' ? similar(A, $elty, ldq, n) : similar(A, $elty, 0)
work = Array($elty, 1)
lwork = BlasInt(-1)
rwork = Array($relty, 2n)
iwork = Array(BlasInt, n)
info = Array(BlasInt, 1)
for i = 1:2
ccall(($(blasfunc(f)), liblapack), Void,
(Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt},
Ptr{BlasInt}),
&jobu, &jobv, &jobq, &m,
&n, &p, k, l,
A, &lda, B, &ldb,
alpha, beta, U, &ldu,
V, &ldv, Q, &ldq,
work, &lwork, rwork, iwork,
info)
@lapackerror
if i == 1
lwork = BlasInt(work[1])
work = Array($elty, lwork)
end
end
if m - k[1] - l[1] >= 0
R = triu(A[1:k[1] + l[1],n - k[1] - l[1] + 1:n])
else
R = triu([A[1:m, n - k[1] - l[1] + 1:n]; B[m - k[1] + 1:l[1], n - k[1] - l[1] + 1:n]])
end
return U, V, Q, alpha, beta, k[1], l[1], R
end
end
end

"""
ggsvd3!(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R)
Finds the generalized singular value decomposition of `A` and `B`, `U'*A*Q = D1*R`
and `V'*B*Q = D2*R`. `D1` has `alpha` on its diagonal and `D2` has `beta` on its
diagonal. If `jobu = U`, the orthogonal/unitary matrix `U` is computed. If
`jobv = V` the orthogonal/unitary matrix `V` is computed. If `jobq = Q`,
the orthogonal/unitary matrix `Q` is computed. If `jobu`, `jobv`, or `jobq` is
`N`, that matrix is not computed. This function requires LAPACK 3.6.0.
"""
ggsvd3!

## Expert driver and generalized eigenvalue problem
for (geevx, ggev, elty) in
((:dgeevx_,:dggev_,:Float64),
Expand All @@ -1679,6 +1837,7 @@ for (geevx, ggev, elty) in
# $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
# $ WI( * ), WORK( * ), WR( * )
function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix{$elty})
checkfinite(A) # balancing routines don't support NaNs and Infs
n = chksquare(A)
lda = max(1,stride(A,2))
wr = similar(A, $elty, n)
Expand Down Expand Up @@ -1825,6 +1984,7 @@ for (geevx, ggev, elty, relty) in
# COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
# $ W( * ), WORK( * )
function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix{$elty})
checkfinite(A) # balancing routines don't support NaNs and Infs
n = chksquare(A)
lda = max(1,stride(A,2))
w = similar(A, $elty, n)
Expand Down Expand Up @@ -4457,6 +4617,7 @@ for (gehrd, elty) in
# * .. Array Arguments ..
# DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
chkstride1(A)
checkfinite(A) # balancing routines don't support NaNs and Infs
n = chksquare(A)
tau = similar(A, $elty, max(0,n - 1))
work = Array($elty, 1)
Expand Down
14 changes: 12 additions & 2 deletions base/linalg/svd.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,12 @@ end
GeneralizedSVD{T}(U::AbstractMatrix{T}, V::AbstractMatrix{T}, Q::AbstractMatrix{T}, a::Vector, b::Vector, k::Int, l::Int, R::AbstractMatrix{T}) = GeneralizedSVD{T,typeof(U)}(U, V, Q, a, b, k, l, R)

function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
U, V, Q, a, b, k, l, R = LAPACK.ggsvd!('U', 'V', 'Q', A, B)
# xggsvd3 replaced xggsvd in LAPACK 3.6.0
if LAPACK.VERSION[] < v"3.6.0"
U, V, Q, a, b, k, l, R = LAPACK.ggsvd!('U', 'V', 'Q', A, B)
else
U, V, Q, a, b, k, l, R = LAPACK.ggsvd3!('U', 'V', 'Q', A, B)
end
GeneralizedSVD(U, V, Q, a, b, Int(k), Int(l), R)
end
svdfact{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = svdfact!(copy(A),copy(B))
Expand Down Expand Up @@ -130,7 +135,12 @@ function getindex{T}(obj::GeneralizedSVD{T}, d::Symbol)
end

function svdvals!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
_, _, _, a, b, k, l, _ = LAPACK.ggsvd!('N', 'N', 'N', A, B)
# xggsvd3 replaced xggsvd in LAPACK 3.6.0
if LAPACK.VERSION[] < v"3.6.0"
_, _, _, a, b, k, l, _ = LAPACK.ggsvd!('N', 'N', 'N', A, B)
else
_, _, _, a, b, k, l, _ = LAPACK.ggsvd3!('N', 'N', 'N', A, B)
end
a[1:k + l] ./ b[1:k + l]
end
svdvals{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = svdvals!(copy(A),copy(B))
Expand Down
8 changes: 7 additions & 1 deletion doc/stdlib/linalg.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1578,7 +1578,13 @@ set of functions in future releases.

.. Docstring generated from Julia source
Finds the generalized singular value decomposition of ``A`` and ``B``\ , ``U'*A*Q = D1*R`` and ``V'*B*Q = D2*R``\ . ``D1`` has ``alpha`` on its diagonal and ``D2`` has ``beta`` on its diagonal. If ``jobu = U``\ , the orthogonal/unitary matrix ``U`` is computed. If ``jobv = V`` the orthogonal/unitary matrix ``V`` is computed. If ``jobq = Q``\ , the orthogonal/unitary matrix ``Q`` is computed. If ``job{u,v,q} = N``\ , that matrix is not computed.
Finds the generalized singular value decomposition of ``A`` and ``B``\ , ``U'*A*Q = D1*R`` and ``V'*B*Q = D2*R``\ . ``D1`` has ``alpha`` on its diagonal and ``D2`` has ``beta`` on its diagonal. If ``jobu = U``\ , the orthogonal/unitary matrix ``U`` is computed. If ``jobv = V`` the orthogonal/unitary matrix ``V`` is computed. If ``jobq = Q``\ , the orthogonal/unitary matrix ``Q`` is computed. If ``jobu``\ , ``jobv`` or ``jobq`` is ``N``\ , that matrix is not computed. This function is only available in LAPACK versions prior to 3.6.0.

.. function:: ggsvd3!(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R)

.. Docstring generated from Julia source
Finds the generalized singular value decomposition of ``A`` and ``B``\ , ``U'*A*Q = D1*R`` and ``V'*B*Q = D2*R``\ . ``D1`` has ``alpha`` on its diagonal and ``D2`` has ``beta`` on its diagonal. If ``jobu = U``\ , the orthogonal/unitary matrix ``U`` is computed. If ``jobv = V`` the orthogonal/unitary matrix ``V`` is computed. If ``jobq = Q``\ , the orthogonal/unitary matrix ``Q`` is computed. If ``jobu``\ , ``jobv``\ , or ``jobq`` is ``N``\ , that matrix is not computed. This function requires LAPACK 3.6.0.

.. function:: geevx!(balanc, jobvl, jobvr, sense, A) -> (A, w, VL, VR, ilo, ihi, scale, abnrm, rconde, rcondv)

Expand Down
17 changes: 16 additions & 1 deletion test/linalg/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,12 @@ for elty in (Float32, Float64, Complex64, Complex128)
@test_approx_eq S lS
@test_approx_eq V' lVt
B = rand(elty,10,10)
@test_throws DimensionMismatch LAPACK.ggsvd!('S','S','S',A,B)
# xggsvd3 replaced xggsvd in LAPACK 3.6.0
if LAPACK.VERSION[] < v"3.6.0"
@test_throws DimensionMismatch LAPACK.ggsvd!('S','S','S',A,B)
else
@test_throws DimensionMismatch LAPACK.ggsvd3!('S','S','S',A,B)
end
end

#geevx, ggev errors
Expand Down Expand Up @@ -492,3 +497,13 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
@test_approx_eq FJulia.τ FLAPACK[2]
end
end

# Issue 13976
let A = [NaN 0.0 NaN; 0 0 0; NaN 0 NaN]
@test_throws ArgumentError expm(A)
end

# Issue 14065 (and 14220)
let A = [NaN NaN; NaN NaN]
@test_throws ArgumentError eigfact(A)
end

2 comments on commit 03c072d

@tkelman
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@nanosoldier runbenchmarks(ALL, vs="@2ac304dfba75fad148d4070ef4f8a2e400c305bb")

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

Please sign in to comment.