Skip to content

Commit

Permalink
Merge pull request #3298 from JuliaLang/anj/gelsy
Browse files Browse the repository at this point in the history 
Wrap xgelsy and change default least square solver from xgelsd to xgelsy
  • Loading branch information
@andreasnoack
andreasnoack committed Jun 5, 2013
2 parents 026f394 + 1d3ac13 commit 9a771d8
Show file tree
Hide file tree
Showing 2 changed files with 107 additions and 21 deletions.
3 changes: 2 additions & 1 deletion base/linalg/dense.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -405,8 +405,9 @@ function (\){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedVecOrMat{T})
ans, _, _, info = LAPACK.gesv!(copy(A), copy(B))
if info > 0; throw(SingularException(info)); end
return ans
else
LAPACK.gelsy!(copy(A), copy(B))[1]
end
LAPACK.gelsd!(copy(A), copy(B))[1]
end

(\){T1<:BlasFloat, T2<:BlasFloat}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) =
Expand Down
125 changes: 105 additions & 20 deletions base/linalg/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -455,8 +455,9 @@ for (gels, gesv, getrs, getri, elty) in
end
end
end
for (gelsd, elty) in ((:dgelsd_, :Float64),
(:sgelsd_, :Float32))
for (gelsd, gelsy, elty) in
((:dgelsd_,:dgelsy_,:Float64),
(:sgelsd_,:sgelsy_,:Float32))
@eval begin
# SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
# $ WORK, LWORK, IWORK, INFO )
Expand All @@ -467,17 +468,13 @@ for (gelsd, elty) in ((:dgelsd_, :Float64),
# * .. Array Arguments ..
# INTEGER IWORK( * )
# DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond)
function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real)
chkstride1(A, B)
m, n = size(A)
if size(B, 1) != m; throw(DimensionMismatch("gelsd!")); end
if size(B, 1) < n
newB = Array($elty, n, size(B, 2))
newB[1:size(B, 1), :] = B
else
newB = B
end
newB = [B; zeros($elty, max(0, n - size(B, 1)), size(B, 2))]
s = Array($elty, min(m, n))
rcond = convert($elty, rcond)
rnk = Array(BlasInt, 1)
info = Array(BlasInt, 1)
work = Array($elty, 1)
Expand All @@ -500,11 +497,56 @@ for (gelsd, elty) in ((:dgelsd_, :Float64),
end
isa(B, Vector) ? newB[1:n] : newB[1:n,:], rnk[1]
end
gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsd!(A, B, -1.)
gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsd!(A, B, -one($elty))

# SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
# $ WORK, LWORK, INFO )
# * .. Scalar Arguments ..
# INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
# DOUBLE PRECISION RCOND
# * ..
# * .. Array Arguments ..
# INTEGER JPVT( * )
# DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
function gelsy!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real)
chkstride1(A, B)
m = size(A, 1)
n = size(A, 2)
nrhs = size(B, 2)
if m != size(B, 1) throw(DimensionMismatch("Left and right hands side must have same number of rows")) end
newB = [B; zeros($elty, max(0, n - size(B, 1)), size(B, 2))]
lda = max(1, m)
ldb = max(1, m, n)
jpvt = Array(BlasInt, n)
rcond = convert($elty, rcond)
rnk = Array(BlasInt, 1)
work = Array($elty, 1)
lwork = -1
info = Array(BlasInt, 1)
for i = 1:2
ccall(($(string(gelsy)), liblapack), Void,
(Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{BlasInt}),
&m, &n, &nrhs, A,
&lda, newB, &ldb, jpvt,
&rcond, rnk, work, &lwork,
info)
if i == 1
lwork = blas_int(work[1])
work = Array($elty, lwork)
end
end
if info[1] != 0 throw(LAPACKException(info[1])) end
return isa(B, Vector) ? newB[1:n] : newB[1:n,:], rnk[1]
end
gelsy!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsy!(A, B, eps($elty))
end
end
for (gelsd, elty, relty) in ((:zgelsd_, :Complex128, :Float64),
(:cgelsd_, :Complex64, :Float32))
for (gelsd, gelsy, elty, relty) in
((:zgelsd_,:zgelsy_,:Complex128,:Float64),
(:cgelsd_,:cgelsy_,:Complex64,:Float32))
@eval begin
# SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
# $ WORK, LWORK, RWORK, IWORK, INFO )
Expand All @@ -516,17 +558,13 @@ for (gelsd, elty, relty) in ((:zgelsd_, :Complex128, :Float64),
# INTEGER IWORK( * )
# DOUBLE PRECISION RWORK( * ), S( * )
# COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond)
function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real)
chkstride1(A, B)
m, n = size(A)
if size(B,1) != m; throw(DimensionMismatch("gelsd!")); end
if size(B, 1) < n
newB = Array($elty, n, size(B, 2))
newB[1:size(B, 1), :] = B
else
newB = B
end
newB = [B; zeros($elty, max(0, n - size(B, 1)), size(B, 2))]
s = Array($elty, min(m, n))
rcond = convert($relty, rcond)
rnk = Array(BlasInt, 1)
info = Array(BlasInt, 1)
work = Array($elty, 1)
Expand All @@ -551,9 +589,56 @@ for (gelsd, elty, relty) in ((:zgelsd_, :Complex128, :Float64),
end
isa(B, Vector) ? newB[1:n] : newB[1:n,:], rnk[1]
end
gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsd!(A, B, -1.)
gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsd!(A, B, -one($relty))

# SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
# $ WORK, LWORK, RWORK, INFO )
# * .. Scalar Arguments ..
# INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
# DOUBLE PRECISION RCOND
# * ..
# * .. Array Arguments ..
# INTEGER JPVT( * )
# DOUBLE PRECISION RWORK( * )
# COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
function gelsy!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real)
chkstride1(A, B)
m = size(A, 1)
n = size(A, 2)
nrhs = size(B, 2)
if m != size(B, 1) throw(DimensionMismatch("Left and right hands side must have same number of rows")) end
newB = [B; zeros($elty, max(0, n - size(B, 1)), size(B, 2))]
lda = max(1, m)
ldb = max(1, m, n)
jpvt = Array(BlasInt, n)
rcond = convert($relty, rcond)
rnk = Array(BlasInt, 1)
work = Array($elty, 1)
lwork = -1
rwork = Array($relty, 2n)
info = Array(BlasInt, 1)
for i = 1:2
ccall(($(string(gelsy)), liblapack), Void,
(Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$relty}, Ptr{BlasInt}),
&m, &n, &nrhs, A,
&lda, newB, &ldb, jpvt,
&rcond, rnk, work, &lwork,
rwork, info)
if i == 1
lwork = blas_int(real(work[1]))
work = Array($elty, lwork)
end
end
if info[1] != 0 throw(LAPACKException(info[1])) end
return isa(B, Vector) ? newB[1:n] : newB[1:n,:], rnk[1]
end
gelsy!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsy!(A, B, eps($relty))
end
end

for (gglse, elty) in ((:dgglse_, :Float64),
(:sgglse_, :Float32),
(:zgglse_, :Complex128),
Expand Down

0 comments on commit 9a771d8

Please sign in to comment.