Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add julia lu factorization. Fix some promotion for linear systems and fix bidiagonal and triangular solvers for multiple right hand side. #5381

Merged
merged 1 commit into from
Jan 17, 2014
Merged

Add julia lu factorization. Fix some promotion for linear systems and fix bidiagonal and triangular solvers for multiple right hand side. #5381

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 21 additions & 14 deletions base/linalg/bidiag.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,36 +117,43 @@ end


#Linear solvers
\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), b::AbstractVector{T}) = naivesub!(A, b, similar(b))
\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...)
function \{TA<:Number,Tb<:Number}(A::Union(Bidiagonal{TA},Triangular{TA}), b::AbstractVector{Tb})
TAb = typeof(one(TA)/one(Tb))
naivesub!(A, b, similar(b, TAb))
end
function \{TA<:Number,TB<:Number}(A::Union(Bidiagonal{TA},Triangular{TA}), B::AbstractMatrix{TB})
TAB = typeof(one(TA)/one(TB))
hcat([naivesub!(A, B[:,i], similar(B[:,i], TAB)) for i=1:size(B,2)]...)
end
A_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(A,b)
At_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(transpose(A),b)
Ac_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(ctranspose(A),b)
for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
function ($func)(A::Bidiagonal, B::AbstractMatrix)
function ($func)(A::Union(Bidiagonal,Triangular), B::AbstractMatrix)
tmp = similar(B[:,1])
n = size(B, 1)
for i=1:size(B,2)
($func)(A, B[:,i])
copy!(tmp, 1, B, (i-1)*n+1, n)
($func)(A, tmp)
copy!(B, (i-1)*n+1, tmp, 1, n) # Modify this when array view are implemented.
end
B
end
end end
for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
function ($func)(A::Bidiagonal, B::AbstractMatrix)
function ($func)(A::Union(Bidiagonal,Triangular), B::AbstractMatrix)
tmp = similar(B[:,2])
m, n = size(B)
nm = n*m
for i=1:size(B,1)
($func)(A, B[i,:])
copy!(tmp, 1, B, i:m:nm, n)
($func)(A, tmp)
copy!(B, i:m:nm, tmp, 1, n)
end
B
end
end end

function inv(A::Union(Bidiagonal, Triangular))
B = eye(eltype(A), size(A, 1))
for i=1:size(B,2)
naivesub!(A, B[:,i])
end
B
end

#Generic solver using naive substitution
function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector=b)
N = size(A, 2)
Expand Down
23 changes: 15 additions & 8 deletions base/linalg/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -412,12 +412,24 @@ function factorize!{T}(A::Matrix{T})
return Triangular(A, :U)
end
if herm
C, info = T <: BlasFloat ? LAPACK.potrf!('U', copy(A)) : LAPACK.potrf!('U', float(A))
if T <: BlasFloat
C, info = LAPACK.potrf!('U', copy(A))
elseif typeof(one(T)/one(T)) <: BlasFloat
C, info = LAPACK.potrf!('U', float(A))
else
error("Unable to factorize hermitian $(typeof(A)). Try converting to other element type or use explicit factorization.")
end
if info == 0 return Cholesky(C, 'U') end
return factorize!(Hermitian(A))
end
if sym
C, info = T <: BlasFloat ? LAPACK.potrf!('U', copy(A)) : LAPACK.potrf!('U', float(A))
if T <: BlasFloat
C, info = LAPACK.potrf!('U', copy(A))
elseif eltype(one(T)/one(T)) <: BlasFloat
C, info = LAPACK.potrf!('U', float(A))
else
error("Unable to factorize symmetric $(typeof(A)). Try converting to other element type or use explicit factorization.")
end
if info == 0 return Cholesky(C, 'U') end
return factorize!(Symmetric(A))
end
Expand All @@ -428,13 +440,8 @@ end

factorize(A::AbstractMatrix) = factorize!(copy(A))

(\){T1<:BlasFloat, T2<:BlasFloat}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) =
(\)(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}) = 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})
function (\)(A::StridedMatrix, B::StridedVecOrMat)
m, n = size(A)
if m == n
if istril(A)
Expand Down
2 changes: 2 additions & 0 deletions base/linalg/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@ type Diagonal{T} <: AbstractMatrix{T}
end
Diagonal(A::Matrix) = Diagonal(diag(A))

convert{T<:Number,S<:Number}(::Type{Diagonal{T}}, A::Diagonal{S}) = T == S ? A : Diagonal(convert(Vector{T}, A.diag))

size(D::Diagonal) = (length(D.diag),length(D.diag))
size(D::Diagonal,d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(D.diag) : 1)

Expand Down
Loading