Skip to content

Commit

Permalink
Merge pull request #5255 from JuliaLang/cjh/solve-triu
Browse files Browse the repository at this point in the history 
More specialized routines for Triangular matrices
  • Loading branch information
@jiahao
jiahao committed Jan 1, 2014
2 parents 26e837b + f15af7a commit 48e599e
Show file tree
Hide file tree
Showing 6 changed files with 473 additions and 27 deletions.
10 changes: 9 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,14 @@ Library improvements

* `sparse(A) \ B` now supports a matrix `B` of right-hand sides ([#5196]).

* More routines for specialized matrix types
- new algorithms for linear solvers and eigensystems of `Triangular`
matrices of generic types ([#5255])
- specialized methods `transpose`, `ctranspose`, `istril`, `istriu` for
`Triangular` ([#5255])
- new LAPACK wrappers
- condition number estimate `cond(A::Triangular)` ([#5255])

Deprecated or removed
---------------------

Expand All @@ -90,7 +98,7 @@ Deprecated or removed
[#4670]: https://github.com/JuliaLang/julia/issues/4670
[#5007]: https://github.com/JuliaLang/julia/issues/5007
[#5076]: https://github.com/JuliaLang/julia/issues/5076

[#5255]: https://github.com/JuliaLang/julia/issues/5255

Julia v0.2.0 Release Notes
==========================
Expand Down
4 changes: 2 additions & 2 deletions base/linalg/dense.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -430,7 +430,7 @@ factorize(A::AbstractMatrix) = factorize!(copy(A))
(\)(convert(Array{promote_type(T1,T2)},A), convert(Array{promote_type(T1,T2)},B))
(\){T1<:BlasFloat, T2<:Number}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(A, convert(Array{T1}, B))
(\){T1<:Number, T2<:BlasFloat}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(convert(Array{T2}, A), B)
(\){T1<:Number, T2<:Number}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(float64(A), float64(B))
(\){T1<:Number, T2<:Number}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = T1<:Complex||T2<:Complex ? (\)(complex128(A), complex128(B)) : (\)(float64(A), float64(B))
(\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B)
function (\){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedVecOrMat{T})
m, n = size(A)
Expand Down Expand Up @@ -473,7 +473,7 @@ function cond(A::StridedMatrix, p::Real=2)
if p == 2
v = svdvals(A)
maxv = maximum(v)
return maxv == 0.0 ? Inf : maxv / minimum(v)
return maxv == 0.0 ? inf(typeof(real(A[1,1]))) : maxv / minimum(v)
elseif p == 1 || p == Inf
chksquare(A)
return cond(lufact(A), p)
Expand Down
2 changes: 1 addition & 1 deletion base/linalg/factorization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,7 @@ Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular

inv(A::LU)=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info

cond(A::LU, p) = 1.0/LinAlg.LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))
cond(A::LU, p) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))

## QR decomposition without column pivots. By the faster geqrt3
type QR{S<:BlasFloat} <: Factorization{S}
Expand Down
230 changes: 228 additions & 2 deletions base/linalg/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1810,8 +1810,8 @@ for (trtri, trtrs, elty) in
chkstride1(A)
n = chksquare(A)
@chkuplo
if size(B,1) != n throw(DimensionMismatch("trtrs!")) end
info = Array(BlasInt, 1)
size(B,1)==n || throw(DimensionMismatch(""))
info = Array(BlasInt, 1)
ccall(($(string(trtrs)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
Expand All @@ -1823,6 +1823,232 @@ for (trtri, trtrs, elty) in
end
end

#Eigenvector computation and condition number estimation
for (trcon, trevc, trrfs, elty) in
((:dtrcon_,:dtrevc_,:dtrrfs_,:Float64),
(:strcon_,:strevc_,:strrfs_,:Float32))
@eval begin
#SUBROUTINE DTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND, WORK,
# IWORK, INFO )
#.. Scalar Arguments ..
#CHARACTER DIAG, NORM, UPLO
#INTEGER INFO, LDA, N
#DOUBLE PRECISION RCOND
#.. Array Arguments ..
#INTEGER IWORK( * )
#DOUBLE PRECISION A( LDA, * ), WORK( * )
function trcon!(norm::BlasChar, uplo::BlasChar, diag::BlasChar,
A::StridedMatrix{$elty})
chkstride1(A)
n = chksquare(A)
@chkuplo
rcond = Array($elty, 1)
work = Array($elty, 3n)
iwork = Array(BlasInt, n)
info = Array(BlasInt, 1)
ccall(($(string(trcon)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
&norm, &uplo, &diag, &n,
A, &max(1,stride(A,2)), rcond, work, iwork, info)
@lapackerror
rcond[1]
end
# SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
# LDVR, MM, M, WORK, INFO )
#
# .. Scalar Arguments ..
# CHARACTER HOWMNY, SIDE
# INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
# ..
# .. Array Arguments ..
# LOGICAL SELECT( * )
# DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
#$ WORK( * )
function trevc!(side::BlasChar, howmny::BlasChar,
select::Vector{Bool}, A::StridedMatrix{$elty},
VL::StridedMatrix{$elty}=similar(A), VR::StridedMatrix{$elty}=similar(A))
chkstride1(A)
chksquare(A)
ldt, n = size(A)
ldvl, mm = size(VL)
ldvr, mm = size(VR)
m = Array(BlasInt, 1)
work = Array($elty, 3n)
info = Array(BlasInt, 1)
ccall(($(string(trevc)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Bool}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
&side, &howmny, select, &n, A, &ldt, VL, &ldvl, VR, &ldvr, &mm,
m, work, info)
@lapackerror

#Decide what exactly to return
if howmny=='S' #compute selected eigenvectors
if side=='L' #left eigenvectors only
return select, VL[:,1:m[1]]
elseif side=='R' #right eigenvectors only
return select, VR[:,1:m[1]]
else #side=='B' #both eigenvectors
return select, VL[:,1:m[1]], VR[:,1:m[1]]
end
else #compute all eigenvectors
if side=='L' #left eigenvectors only
return VL[:,1:m[1]]
elseif side=='R' #right eigenvectors only
return VR[:,1:m[1]]
else #side=='B' #both eigenvectors
return VL[:,1:m[1]], VR[:,1:m[1]]
end
end
end
# SUBROUTINE DTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
# LDX, FERR, BERR, WORK, IWORK, INFO )
# .. Scalar Arguments ..
# CHARACTER DIAG, TRANS, UPLO
# INTEGER INFO, LDA, LDB, LDX, N, NRHS
# .. Array Arguments ..
# INTEGER IWORK( * )
# DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
#$ WORK( * ), X( LDX, * )
function trrfs!(uplo::BlasChar, trans::BlasChar, diag::BlasChar,
A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, X::StridedVecOrMat{$elty},
Ferr::StridedVector{$elty}=Array($elty, size(B,2)), Berr::StridedVector{$elty}=Array($elty, size(B,2)))
@chkuplo
n=size(A,2)
nrhs=size(B,2)
nrhs==size(X,2) || throw(DimensionMismatch(""))
work=Array($elty, 3n)
iwork=Array(BlasInt, n)
info=Array(BlasInt, 1)
ccall(($(string(trrfs)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
&uplo, &trans, &diag, &n,
&nrhs, A, &max(1,stride(A,2)), B, &max(1,stride(B,2)), X, &max(1,stride(X,2)),
Ferr, Berr, work, iwork, info)
@lapackerror
Ferr, Berr
end
end
end
for (trcon, trevc, trrfs, elty, relty) in
((:ztrcon_,:ztrevc_,:ztrrfs_,:Complex128,:Float64),
(:ctrcon_,:ctrevc_,:ctrrfs_,:Complex64, :Float32))
@eval begin
#SUBROUTINE ZTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND, WORK,
# RWORK, INFO )
#.. Scalar Arguments ..
#CHARACTER DIAG, NORM, UPLO
#INTEGER INFO, LDA, N
#DOUBLE PRECISION RCOND
#.. Array Arguments ..
#DOUBLE PRECISION RWORK( * )
#COMPLEX*16 A( LDA, * ), WORK( * )
function trcon!(norm::BlasChar, uplo::BlasChar, diag::BlasChar,
A::StridedMatrix{$elty})
chkstride1(A)
n = chksquare(A)
@chkuplo
rcond = Array($relty, 1)
work = Array($elty, 2n)
rwork = Array($elty, n)
info = Array(BlasInt, 1)
ccall(($(string(trcon)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}),
&norm, &uplo, &diag, &n,
A, &max(1,stride(A,2)), rcond, work, rwork, info)
@lapackerror
rcond[1]
end

# SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
# LDVR, MM, M, WORK, RWORK, INFO )
#
# .. Scalar Arguments ..
# CHARACTER HOWMNY, SIDE
# INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
# ..
# .. Array Arguments ..
# LOGICAL SELECT( * )
# DOUBLE PRECISION RWORK( * )
# COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
#$ WORK( * )
function trevc!(side::BlasChar, howmny::BlasChar,
select::Vector{Bool}, A::StridedMatrix{$elty},
VL::StridedMatrix{$elty}=similar(A), VR::StridedMatrix{$elty}=similar(A))
chkstride1(A)
chksquare(A)
ldt, n = size(A)
ldvl, mm = size(VL)
ldvr, mm = size(VR)
m = Array(BlasInt, 1)
work = Array($elty, 2n)
rwork = Array($relty, n)
info = Array(BlasInt, 1)
ccall(($(string(trevc)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Bool}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}),
&side, &howmny, select, &n, A, &ldt, VL, &ldvl, VR, &ldvr, &mm,
m, work, rwork, info)
@lapackerror

#Decide what exactly to return
if howmny=='S' #compute selected eigenvectors
if side=='L' #left eigenvectors only
return select, VL[:,1:m[1]]
elseif side=='R' #right eigenvectors only
return select, VR[:,1:m[1]]
else #side=='B' #both eigenvectors
return select, VL[:,1:m[1]], VR[:,1:m[1]]
end
else #compute all eigenvectors
if side=='L' #left eigenvectors only
return VL[:,1:m[1]]
elseif side=='R' #right eigenvectors only
return VR[:,1:m[1]]
else #side=='B' #both eigenvectors
return VL[:,1:m[1]], VR[:,1:m[1]]
end
end
end

# SUBROUTINE ZTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
# LDX, FERR, BERR, WORK, IWORK, INFO )
# .. Scalar Arguments ..
# CHARACTER DIAG, TRANS, UPLO
# INTEGER INFO, LDA, LDB, LDX, N, NRHS
# .. Array Arguments ..
# INTEGER IWORK( * )
# DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
#$ WORK( * ), X( LDX, * )
function trrfs!(uplo::BlasChar, trans::BlasChar, diag::BlasChar,
A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, X::StridedVecOrMat{$elty},
Ferr::StridedVector{$relty}=Array($relty, size(B,2)), Berr::StridedVector{$relty}=Array($relty, size(B,2)))
@chkuplo
n=size(A,2)
nrhs=size(B,2)
nrhs==size(X,2) || throw(DimensionMismatch(""))
work=Array($elty, 2n)
rwork=Array($elty, n)
info=Array(BlasInt, 1)
ccall(($(string(trrfs)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}),
&uplo, &trans, &diag, &n,
&nrhs, A, &max(1,stride(A,2)), B, &max(1,stride(B,2)), X, &max(1,stride(X,2)),
Ferr, Berr, work, rwork, info)
@lapackerror
Ferr, Berr
end
end
end

## (ST) Symmetric tridiagonal - eigendecomposition
for (stev, stebz, stegr, stein, elty) in
((:dstev_,:dstebz_,:dstegr_,:dstein_,:Float64),
Expand Down
Loading

0 comments on commit 48e599e

Please sign in to comment.