From e642034d801f254aa8dd4ec7a8e7f02c566cbedf Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 14 Jan 2015 20:39:18 -0500 Subject: [PATCH] Redesign of triangular matrix types. Avoid Upper/lower and unit diagonal specification by parameters. This fixes #9191, type instability of \ for BlasFloat types when no promotion is necessary. --- NEWS.md | 2 + base/exports.jl | 3 +- base/linalg.jl | 3 +- base/linalg/bidiag.jl | 13 +- base/linalg/cholesky.jl | 49 ++- base/linalg/dense.jl | 16 +- base/linalg/diagonal.jl | 3 +- base/linalg/factorization.jl | 8 +- base/linalg/generic.jl | 6 +- base/linalg/lu.jl | 4 +- base/linalg/sparse.jl | 2 +- base/linalg/special.jl | 26 +- base/linalg/triangular.jl | 736 +++++++++++++++++++++++++--------- base/linalg/uniformscaling.jl | 2 +- doc/manual/linear-algebra.rst | 11 +- doc/stdlib/linalg.rst | 2 +- test/linalg/triangular.jl | 278 ++++++------- test/linalg2.jl | 2 +- test/linalg4.jl | 14 +- 19 files changed, 787 insertions(+), 393 deletions(-) diff --git a/NEWS.md b/NEWS.md index 967e03a48513d..3e4a15ceaf01b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -154,6 +154,8 @@ Library improvements * copy(DArray) will now make a copy of the DArray ([#9745]) + * Split `Triangular` type into `UpperTriangular`, `LowerTriangular`, `UnitUpperTriagular` and `UnitLowerTriangular` ([#9779]) + Deprecated or removed --------------------- diff --git a/base/exports.jl b/base/exports.jl index e319785717794..3b257f3e38a94 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -62,6 +62,7 @@ export IOBuffer, IOStream, LocalProcess, + LowerTriangular, MathConst, Matrix, MergeSort, @@ -108,9 +109,9 @@ export SymTridiagonal, Timer, TmStruct, - Triangular, Tridiagonal, UnitRange, + UpperTriangular, UTF16String, UTF32String, Val, diff --git a/base/linalg.jl b/base/linalg.jl index cfc7c74eca088..e81942bb55c1c 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -30,7 +30,8 @@ export SVD, Hermitian, Symmetric, - Triangular, + LowerTriangular, + UpperTriangular, Diagonal, UniformScaling, diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl index b804a5f03a8fe..8701c11b8bdf5 100644 --- a/base/linalg/bidiag.jl +++ b/base/linalg/bidiag.jl @@ -114,24 +114,23 @@ end /(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper) ==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper) -SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular) +SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, TriangularUnion) *(A::SpecialMatrix, B::SpecialMatrix)=full(A)*full(B) #Generic multiplication for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc) @eval begin ($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(full(A), B) - #($func){T}(A::AbstractArray{T}, B::Triangular{T}) = ($func)(full(A), B) end end #Linear solvers -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) +A_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(A, b) +At_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(transpose(A), b) +Ac_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(ctranspose(A), b) for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin - function ($func)(A::Union(Bidiagonal, Triangular), B::AbstractMatrix) + function ($func)(A::Union(Bidiagonal, TriangularUnion), B::AbstractMatrix) tmp = similar(B[:,1]) n = size(B, 1) for i = 1:size(B,2) @@ -143,7 +142,7 @@ for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin end end end for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin - function ($func)(A::Union(Bidiagonal, Triangular), B::AbstractMatrix) + function ($func)(A::Union(Bidiagonal, TriangularUnion), B::AbstractMatrix) tmp = similar(B[:, 2]) m, n = size(B) nm = n*m diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index 4ee7e30e17c4e..1fca762bf4610 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -18,12 +18,35 @@ immutable CholeskyPivoted{T,S<:AbstractMatrix} <: Factorization{T} end CholeskyPivoted{T}(UL::AbstractMatrix{T}, uplo::Char, piv::Vector{BlasInt}, rank::BlasInt, tol::Real, info::BlasInt) = CholeskyPivoted{T,typeof(UL)}(UL, uplo, piv, rank, tol, info) -function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U) +function chol!{T<:BlasFloat}(A::StridedMatrix{T}) + C, info = LAPACK.potrf!('U', A) + return @assertposdef UpperTriangular(C) info +end +function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol) C, info = LAPACK.potrf!(char_uplo(uplo), A) - return @assertposdef Triangular{eltype(C),typeof(C),uplo,false}(C) info + return @assertposdef uplo == :U ? UpperTriangular(C) : LowerTriangular(C) info end -function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U) +function chol!{T}(A::AbstractMatrix{T}) + n = chksquare(A) + @inbounds begin + for k = 1:n + for i = 1:k - 1 + A[k,k] -= A[i,k]'A[i,k] + end + A[k,k] = chol!(A[k,k], uplo) + AkkInv = inv(A[k,k]) + for j = k + 1:n + for i = 1:k - 1 + A[k,j] -= A[i,k]'A[i,j] + end + A[k,j] = A[k,k]'\A[k,j] + end + end + end + return UpperTriangular(A) +end +function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol) n = chksquare(A) @inbounds begin if uplo == :L @@ -58,7 +81,7 @@ function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U) throw(ArgumentError("uplo must be either :U or :L but was $(uplo)")) end end - return Triangular(A, uplo, false) + return uplo == :U ? UpperTriangular(A) : LowerTriangular(A) end function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) @@ -118,14 +141,14 @@ size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL) size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d) function getindex{T,S,UpLo}(C::Cholesky{T,S,UpLo}, d::Symbol) - d == :U && return Triangular(UpLo == d ? C.UL : C.UL',:U) - d == :L && return Triangular(UpLo == d ? C.UL : C.UL',:L) - d == :UL && return Triangular(C.UL, UpLo) + d == :U && return UpperTriangular(UpLo == d ? C.UL : C.UL') + d == :L && return LowerTriangular(UpLo == d ? C.UL : C.UL') + d == :UL && return UpLo == :U ? UpperTriangular(C.UL) : LowerTriangular(C.UL) throw(KeyError(d)) end function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol) - d == :U && return Triangular(symbol(C.uplo) == d ? C.UL : C.UL', :U) - d == :L && return Triangular(symbol(C.uplo) == d ? C.UL : C.UL', :L) + d == :U && return UpperTriangular(symbol(C.uplo) == d ? C.UL : C.UL') + d == :L && return LowerTriangular(symbol(C.uplo) == d ? C.UL : C.UL') d == :p && return C.piv if d == :P n = size(C, 1) @@ -142,8 +165,8 @@ show{T,S<:AbstractMatrix,UpLo}(io::IO, C::Cholesky{T,S,UpLo}) = (println("$(type A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat{T}) = LAPACK.potrs!('U', C.UL, B) A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat{T}) = LAPACK.potrs!('L', C.UL, B) -A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat) = Ac_ldiv_B!(Triangular(C.UL, :L, false), A_ldiv_B!(Triangular(C.UL, :L, false), B)) -A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat) = A_ldiv_B!(Triangular(C.UL, :U, false), Ac_ldiv_B!(Triangular(C.UL, :U, false), B)) +A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat) = Ac_ldiv_B!(LowerTriangular(C.UL), A_ldiv_B!(LowerTriangular(C.UL), B)) +A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat) = A_ldiv_B!(UpperTriangular(C.UL), Ac_ldiv_B!(UpperTriangular(C.UL), B)) function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T}) chkfullrank(C) @@ -161,8 +184,8 @@ function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T}) end B end -A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv]))[invperm(C.piv)] -A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv,:]))[invperm(C.piv),:] +A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(LowerTriangular(C.UL), A_ldiv_B!(LowerTriangular(C.UL), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(UpperTriangular(C.UL), Ac_ldiv_B!(UpperTriangular(C.UL), B[C.piv]))[invperm(C.piv)] +A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(LowerTriangular(C.UL), A_ldiv_B!(LowerTriangular(C.UL), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(UpperTriangular(C.UL), Ac_ldiv_B!(UpperTriangular(C.UL), B[C.piv,:]))[invperm(C.piv),:] function det{T,S,UpLo}(C::Cholesky{T,S,UpLo}) dd = one(T) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index 2a1966753c5e4..5477c91a03afb 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -307,7 +307,7 @@ function sqrtm{T<:Real}(A::StridedMatrix{T}) issym(A) && return sqrtm(Symmetric(A)) n = chksquare(A) SchurF = schurfact(complex(A)) - R = full(sqrtm(Triangular(SchurF[:T], :U, false))) + R = full(sqrtm(UpperTriangular(SchurF[:T]))) retmat = SchurF[:vectors]*R*SchurF[:vectors]' all(imag(retmat) .== 0) ? real(retmat) : retmat end @@ -315,7 +315,7 @@ function sqrtm{T<:Complex}(A::StridedMatrix{T}) ishermitian(A) && return sqrtm(Hermitian(A)) n = chksquare(A) SchurF = schurfact(A) - R = full(sqrtm(Triangular(SchurF[:T], :U, false))) + R = full(sqrtm(UpperTriangular(SchurF[:T]))) SchurF[:vectors]*R*SchurF[:vectors]' end sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b) @@ -325,9 +325,9 @@ function inv{S}(A::StridedMatrix{S}) T = typeof(one(S)/one(S)) Ac = convert(AbstractMatrix{T}, A) if istriu(Ac) - Ai = inv(Triangular(A, :U, false)) + Ai = inv(UpperTriangular(A)) elseif istril(Ac) - Ai = inv(Triangular(A, :L, false)) + Ai = inv(LowerTriangular(A)) else Ai = inv(lufact(Ac)) end @@ -377,7 +377,7 @@ function factorize{T}(A::Matrix{T}) if utri1 return Bidiagonal(diag(A), diag(A, -1), false) end - return Triangular(A, :L) + return LowerTriangular(A) end if utri return Bidiagonal(diag(A), diag(A, 1), true) @@ -392,7 +392,7 @@ function factorize{T}(A::Matrix{T}) end end if utri - return Triangular(A, :U) + return UpperTriangular(A) end if herm try @@ -413,9 +413,9 @@ function (\)(A::StridedMatrix, B::StridedVecOrMat) m, n = size(A) if m == n if istril(A) - return istriu(A) ? \(Diagonal(A),B) : \(Triangular(A, :L),B) + return istriu(A) ? \(Diagonal(A),B) : \(LowerTriangular(A),B) end - istriu(A) && return \(Triangular(A, :U),B) + istriu(A) && return \(UpperTriangular(A),B) return \(lufact(A),B) end return qrfact(A,pivot=eltype(A)<:BlasFloat)\B diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl index 7a621e8a73ec5..d2c52e328271b 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/diagonal.jl @@ -8,7 +8,8 @@ Diagonal(A::Matrix) = Diagonal(diag(A)) convert{T}(::Type{Diagonal{T}}, D::Diagonal{T}) = D convert{T}(::Type{Diagonal{T}}, D::Diagonal) = Diagonal{T}(convert(Vector{T}, D.diag)) convert{T}(::Type{AbstractMatrix{T}}, D::Diagonal) = convert(Diagonal{T}, D) -convert{T}(::Type{Triangular}, A::Diagonal{T}) = Triangular{T, Diagonal{T}, :U, false}(A) +convert{T}(::Type{UpperTriangular}, A::Diagonal{T}) = UpperTriangular(A) +convert{T}(::Type{LowerTriangular}, A::Diagonal{T}) = LowerTriangular(A) function similar{T}(D::Diagonal, ::Type{T}, d::(Int,Int)) d[1] == d[2] || throw(ArgumentError("Diagonal matrix must be square")) diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 63fe761cc4fa2..9fc00ab461dae 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -274,7 +274,7 @@ function A_mul_Bc!{T}(A::AbstractMatrix{T},Q::QRPackedQ{T}) end A end -A_mul_Bc(A::Triangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) +A_mul_Bc(A::TriangularUnion, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) function A_mul_Bc{TA,TB}(A::AbstractArray{TA}, B::Union(QRCompactWYQ{TB},QRPackedQ{TB})) TAB = promote_type(TA,TB) A_mul_Bc!(size(A,2)==size(B.factors,1) ? (TA == TAB ? copy(A) : convert(AbstractMatrix{TAB}, A)) : @@ -284,8 +284,8 @@ function A_mul_Bc{TA,TB}(A::AbstractArray{TA}, B::Union(QRCompactWYQ{TB},QRPacke convert(AbstractMatrix{TAB}, B)) end -A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, b::StridedVector{T}) = (A_ldiv_B!(Triangular(A[:R], :U), sub(Ac_mul_B!(A[:Q], b), 1:size(A, 2))); b) -A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, B::StridedMatrix{T}) = (A_ldiv_B!(Triangular(A[:R], :U), sub(Ac_mul_B!(A[:Q], B), 1:size(A, 2), 1:size(B, 2))); B) +A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, b::StridedVector{T}) = (A_ldiv_B!(UpperTriangular(A[:R]), sub(Ac_mul_B!(A[:Q], b), 1:size(A, 2))); b) +A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, B::StridedMatrix{T}) = (A_ldiv_B!(UpperTriangular(A[:R]), sub(Ac_mul_B!(A[:Q], B), 1:size(A, 2), 1:size(B, 2))); B) # Julia implementation similarly to xgelsy function A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) @@ -313,7 +313,7 @@ function A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Re # if cond(r[1:rnk, 1:rnk])*rcond < 1 break end end C, τ = LAPACK.tzrzf!(A.factors[1:rnk,:]) - A_ldiv_B!(Triangular(C[1:rnk,1:rnk],:U),sub(Ac_mul_B!(getq(A),sub(B, 1:mA, 1:nrhs)),1:rnk,1:nrhs)) + A_ldiv_B!(UpperTriangular(C[1:rnk,1:rnk]),sub(Ac_mul_B!(getq(A),sub(B, 1:mA, 1:nrhs)),1:rnk,1:nrhs)) B[rnk+1:end,:] = zero(T) LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, τ, sub(B,1:nA,1:nrhs)) B[1:nA,:] = sub(B, 1:nA, :)[invperm(A[:p]::Vector{BlasInt}),:] diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index 419f947cc2dd6..861d101afe83b 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -1,5 +1,9 @@ ## linalg.jl: Some generic Linear Algebra definitions +# Fall back arithmetic ++(A::AbstractMatrix, B::AbstractMatrix) = full(A) + full(B) +-(A::AbstractMatrix, B::AbstractMatrix) = full(A) - full(B) + scale(X::AbstractArray, s::Number) = scale!(copy(X), s) scale(s::Number, X::AbstractArray) = scale!(copy(X), s) @@ -424,7 +428,7 @@ function elementaryRightTrapezoid!(A::AbstractMatrix, row::Integer) end function det(A::AbstractMatrix) - (istriu(A) || istril(A)) && return det(Triangular(A, :U, false)) + (istriu(A) || istril(A)) && return det(UpperTriangular(A)) return det(lufact(A)) end det(x::Number) = x diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl index 38450cb86b3b7..4fd661735eabc 100644 --- a/base/linalg/lu.jl +++ b/base/linalg/lu.jl @@ -114,8 +114,8 @@ function getindex{T,S<:StridedMatrix}(A::LU{T,S}, d::Symbol) end A_ldiv_B!{T<:BlasFloat, S<:StridedMatrix}(A::LU{T, S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info -A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, b::StridedVector) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), b[ipiv2perm(A.ipiv, length(b))])) -A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, B::StridedMatrix) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), B[ipiv2perm(A.ipiv, size(B, 1)),:])) +A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, b::StridedVector) = A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), b[ipiv2perm(A.ipiv, length(b))])) +A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, B::StridedMatrix) = A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), B[ipiv2perm(A.ipiv, size(B, 1)),:])) At_ldiv_B{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, copy(B)) A.info Ac_ldiv_B{T<:BlasComplex,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, copy(B)) A.info At_ldiv_Bt{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info diff --git a/base/linalg/sparse.jl b/base/linalg/sparse.jl index 14c4c3081c923..d507bbb8e6a53 100644 --- a/base/linalg/sparse.jl +++ b/base/linalg/sparse.jl @@ -126,7 +126,7 @@ end *{TvA,TiA}(X::BitArray{2}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (AbstractMatrix, SparseMatrixCSC), X, A) # TODO: Tridiagonal * Sparse should be implemented more efficiently *{TX,TvA,TiA}(X::Tridiagonal{TX}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (Tridiagonal, AbstractMatrix), X, A) -*{TvA,TiA}(X::Triangular, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A +*{TvA,TiA}(X::TriangularUnion, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A function *{TX,TvA,TiA}(X::AbstractMatrix{TX}, A::SparseMatrixCSC{TvA,TiA}) mX, nX = size(X) nX == A.m || throw(DimensionMismatch()) diff --git a/base/linalg/special.jl b/base/linalg/special.jl index 4113e49446d95..0c0af842fe4c0 100644 --- a/base/linalg/special.jl +++ b/base/linalg/special.jl @@ -4,8 +4,12 @@ convert{T}(::Type{Bidiagonal}, A::Diagonal{T})=Bidiagonal(A.diag, zeros(T, size(A.diag,1)-1), true) convert{T}(::Type{SymTridiagonal}, A::Diagonal{T})=SymTridiagonal(A.diag, zeros(T, size(A.diag,1)-1)) convert{T}(::Type{Tridiagonal}, A::Diagonal{T})=Tridiagonal(zeros(T, size(A.diag,1)-1), A.diag, zeros(T, size(A.diag,1)-1)) -convert(::Type{Triangular}, A::Diagonal) = Triangular(full(A), :L) -convert(::Type{Triangular}, A::Bidiagonal) = Triangular(full(A), A.isupper ? :U : :L) +convert(::Type{UpperTriangular}, A::Diagonal) = UpperTriangular(full(A), :L) +convert(::Type{UnitUpperTriangular}, A::Diagonal) = UnitUpperTriangular(full(A), :L) +convert(::Type{LowerTriangular}, A::Diagonal) = LowerTriangular(full(A), :L) +convert(::Type{UnitLowerTriangular}, A::Diagonal) = UnitLowerTriangular(full(A), :L) +convert(::Type{LowerTriangular}, A::Bidiagonal) = !A.isupper ? LowerTriangular(full(A)) : throw(ArgumentError("Bidiagonal matrix must have lower off diagonal to be converted to LowerTriangular")) +convert(::Type{UpperTriangular}, A::Bidiagonal) = A.isupper ? UpperTriangular(full(A)) : throw(ArgumentError("Bidiagonal matrix must have upper off diagonal to be converted to UpperTriangular")) convert(::Type{Matrix}, D::Diagonal) = diagm(D.diag) function convert(::Type{Diagonal}, A::Union(Bidiagonal, SymTridiagonal)) @@ -42,12 +46,12 @@ function convert(::Type{SymTridiagonal}, A::Tridiagonal) SymTridiagonal(A.d, A.dl) end -function convert(::Type{Diagonal}, A::Triangular) +function convert(::Type{Diagonal}, A::TriangularUnion) full(A) == diagm(diag(A)) || throw(ArgumentError("Matrix cannot be represented as Diagonal")) Diagonal(diag(A)) end -function convert(::Type{Bidiagonal}, A::Triangular) +function convert(::Type{Bidiagonal}, A::TriangularUnion) fA = full(A) if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) return Bidiagonal(diag(A), diag(fA,1), true) @@ -58,9 +62,9 @@ function convert(::Type{Bidiagonal}, A::Triangular) end end -convert(::Type{SymTridiagonal}, A::Triangular) = convert(SymTridiagonal, convert(Tridiagonal, A)) +convert(::Type{SymTridiagonal}, A::TriangularUnion) = convert(SymTridiagonal, convert(Tridiagonal, A)) -function convert(::Type{Tridiagonal}, A::Triangular) +function convert(::Type{Tridiagonal}, A::TriangularUnion) fA = full(A) if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) + diagm(diag(fA, -1), -1) return Tridiagonal(diag(fA, -1), diag(A), diag(fA,1)) @@ -82,7 +86,7 @@ macro commutative(myexpr) end for op in (:+, :-) - SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Triangular, :Matrix] + SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Matrix] for (idx, matrixtype1) in enumerate(SpecialMatrices) #matrixtype1 is the sparser matrix type for matrixtype2 in SpecialMatrices[idx+1:end] #matrixtype2 is the denser matrix type @eval begin #TODO quite a few of these conversions are NOT defined... @@ -93,7 +97,7 @@ for op in (:+, :-) end for matrixtype1 in (:SymTridiagonal,) #matrixtype1 is the sparser matrix type - for matrixtype2 in (:Tridiagonal, :Triangular, :Matrix) #matrixtype2 is the denser matrix type + for matrixtype2 in (:Tridiagonal, :Matrix) #matrixtype2 is the denser matrix type @eval begin ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B) ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B)) @@ -111,7 +115,7 @@ for op in (:+, :-) end end -A_mul_Bc!(A::Triangular, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B) -A_mul_Bc!(A::Triangular, B::QRPackedQ) = A_mul_Bc!(full!(A),B) -A_mul_Bc(A::Triangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) +A_mul_Bc!(A::TriangularUnion, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B) +A_mul_Bc!(A::TriangularUnion, B::QRPackedQ) = A_mul_Bc!(full!(A),B) +A_mul_Bc(A::TriangularUnion, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index b344573187ab1..b0658b1029e3f 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -1,141 +1,226 @@ ## Triangular -immutable Triangular{T,S<:AbstractMatrix,UpLo,IsUnit} <: AbstractMatrix{T} - data::S -end -function Triangular{T}(A::AbstractMatrix{T}, uplo::Symbol, isunit::Bool=false) - chksquare(A) - Triangular{T,typeof(A),uplo,isunit}(A) + +# First loop through all methods that don't need special care for upper/lower and unit diagonal +for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular) + @eval begin + immutable $t{T,S<:AbstractMatrix} <: AbstractMatrix{T} + data::S + end + function $t(A::AbstractMatrix) + Base.LinAlg.chksquare(A) + return $t{eltype(A), typeof(A)}(A) + end + + size(A::$t, args...) = size(A.data, args...) + + convert{T,S}(::Type{$t{T}}, A::$t{T,S}) = A + convert{Tnew,Told,S}(::Type{$t{Tnew}}, A::$t{Told,S}) = (Anew = convert(AbstractMatrix{Tnew}, A.data); $t(Anew)) + convert{Tnew,Told,S}(::Type{AbstractMatrix{Tnew}}, A::$t{Told,S}) = convert($t{Tnew}, A) + convert{T,S}(::Type{Matrix}, A::$t{T,S}) = convert(Matrix{T}, A) + + function similar{T,S,Tnew}(A::$t{T,S}, ::Type{Tnew}, dims::Dims) + dims[1] == dims[2] || throw(ArgumentError("a Triangular matrix must be square")) + length(dims) == 2 || throw(ArgumentError("a Triangular matrix must have two dimensions")) + B = similar(A.data, Tnew, dims) + return $t(B) + end + + copy{T,S}(A::$t{T,S}) = $t{T,S}(copy(A.data)) + + big(A::$t) = $t(big(A.data)) + + real{T<:Real}(A::$t{T}) = A + real{T<:Complex}(A::$t{T}) = (B = real(A.data); $t(B)) + + end end -size(A::Triangular, args...) = size(A.data, args...) +# some cases can be handled with a union +typealias TriangularUnion{T,S} Union(UpperTriangular{T,S}, UnitUpperTriangular{T,S}, LowerTriangular{T,S}, UnitLowerTriangular{T,S}) + +full{T,S}(A::TriangularUnion{T,S}) = convert(Matrix, A) -convert{T,S,UpLo,IsUnit}(::Type{Triangular{T}}, A::Triangular{T,S,UpLo,IsUnit}) = A -convert{Tnew,Told,S,UpLo,IsUnit}(::Type{Triangular{Tnew}}, A::Triangular{Told,S,UpLo,IsUnit}) = (Anew = convert(AbstractMatrix{Tnew}, A.data); Triangular(Anew, UpLo, IsUnit)) -convert{Tnew,Told,S,UpLo,IsUnit}(::Type{AbstractMatrix{Tnew}}, A::Triangular{Told,S,UpLo,IsUnit}) = convert(Triangular{Tnew}, A) -function convert{Tret,T,S,UpLo,IsUnit}(::Type{Matrix{Tret}}, A::Triangular{T,S,UpLo,IsUnit}) +fill!(A::TriangularUnion, x) = (fill!(A.data, x); A) + +# then handle all methods that requires specific handling of upper/lower and unit diagonal + +function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::LowerTriangular{T,S}) B = Array(Tret, size(A, 1), size(A, 1)) copy!(B, A.data) - (UpLo == :L ? tril! : triu!)(B) - if IsUnit - for i = 1:size(B,1) - B[i,i] = 1 - end + tril!(B) + B +end +function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::UnitLowerTriangular{T,S}) + B = Array(Tret, size(A, 1), size(A, 1)) + copy!(B, A.data) + tril!(B) + for i = 1:size(B,1) + B[i,i] = 1 + end + B +end +function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::UpperTriangular{T,S}) + B = Array(Tret, size(A, 1), size(A, 1)) + copy!(B, A.data) + triu!(B) + B +end +function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::UnitUpperTriangular{T,S}) + B = Array(Tret, size(A, 1), size(A, 1)) + copy!(B, A.data) + triu!(B) + for i = 1:size(B,1) + B[i,i] = 1 end B end -convert{T,S,UpLo,IsUnit}(::Type{Matrix}, A::Triangular{T,S,UpLo,IsUnit}) = convert(Matrix{T}, A) -function full!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) +function full!{T,S}(A::LowerTriangular{T,S}) B = A.data - (UpLo == :L ? tril! : triu!)(B) - if IsUnit - for i = 1:size(A,1) - B[i,i] = 1 - end + tril!(B) + B +end +function full!{T,S}(A::UnitLowerTriangular{T,S}) + B = A.data + tril!(B) + for i = 1:size(A,1) + B[i,i] = 1 end B end -full{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = convert(Matrix, A) - -fill!(A::Triangular, x) = (fill!(A.data, x); A) - -function similar{T,S,UpLo,IsUnit,Tnew}(A::Triangular{T,S,UpLo,IsUnit}, ::Type{Tnew}, dims::Dims) - dims[1] == dims[2] || throw(ArgumentError("a Triangular matrix must be square")) - length(dims) == 2 || throw(ArgumentError("a Triangular matrix must have two dimensions")) - A = similar(A.data, Tnew, dims) - return Triangular{Tnew, typeof(A), UpLo, IsUnit}(A) +function full!{T,S}(A::UpperTriangular{T,S}) + B = A.data + triu!(B) + B end - -copy{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T,S,UpLo,IsUnit}(copy(A.data)) - -getindex(A::Triangular, i::Integer) = ((m, n) = divrem(i - 1, size(A, 1)); A[n + 1, m + 1]) -getindex{T,S}(A::Triangular{T,S,:L,true}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(A.data[i,j])) -getindex{T,S}(A::Triangular{T,S,:L,false}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(A.data[i,j]) -getindex{T,S}(A::Triangular{T,S,:U,true}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(A.data[i,j])) -getindex{T,S}(A::Triangular{T,S,:U,false}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(A.data[i,j]) - -istril{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :L -istriu{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :U - -transpose{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(transpose(A.data)) -ctranspose{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(ctranspose(A.data)) -transpose!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(copytri!(A.data, UpLo == :L ? 'L' : 'U')) -ctranspose!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(copytri!(A.data, UpLo == :L ? 'L' : 'U', true)) -diag{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = IsUnit ? ones(T, size(A,1)) : diag(A.data) -function big{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) - M = big(A.data) - Triangular{eltype(M),typeof(M),UpLo,IsUnit}(M) +function full!{T,S}(A::UnitUpperTriangular{T,S}) + B = A.data + triu!(B) + for i = 1:size(A,1) + B[i,i] = 1 + end + B end -real{T<:Real}(A::Triangular{T}) = A -real{T<:Complex,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = (B = real(A.data); Triangular{eltype(B), typeof(B), UpLo, IsUnit}(B)) +getindex(A::TriangularUnion, i::Integer) = ((m, n) = divrem(i - 1, size(A, 1)); A[n + 1, m + 1]) +getindex{T,S}(A::UnitLowerTriangular{T,S}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(A.data[i,j])) +getindex{T,S}(A::LowerTriangular{T,S}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(A.data[i,j]) +getindex{T,S}(A::UnitUpperTriangular{T,S}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(A.data[i,j])) +getindex{T,S}(A::UpperTriangular{T,S}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(A.data[i,j]) + +istril{T,S}(A::LowerTriangular{T,S}) = true +istril{T,S}(A::UnitLowerTriangular{T,S}) = true +istril{T,S}(A::UpperTriangular{T,S}) = false +istril{T,S}(A::UnitUpperTriangular{T,S}) = false +istriu{T,S}(A::LowerTriangular{T,S}) = false +istriu{T,S}(A::UnitLowerTriangular{T,S}) = false +istriu{T,S}(A::UpperTriangular{T,S}) = true +istriu{T,S}(A::UnitUpperTriangular{T,S}) = true + +transpose{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(transpose(A.data)) +transpose{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(transpose(A.data)) +transpose{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(transpose(A.data)) +transpose{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(transpose(A.data)) +ctranspose{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(ctranspose(A.data)) +ctranspose{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(ctranspose(A.data)) +ctranspose{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(ctranspose(A.data)) +ctranspose{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(ctranspose(A.data)) + +transpose!{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(copytri!(A.data, 'L')) +transpose!{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(copytri!(A.data, 'L')) +transpose!{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(copytri!(A.data, 'U')) +transpose!{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(copytri!(A.data, 'U')) +ctranspose!{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(copytri!(A.data, 'L' , true)) +ctranspose!{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(copytri!(A.data, 'L' , true)) +ctranspose!{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(copytri!(A.data, 'U' , true)) +ctranspose!{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(copytri!(A.data, 'U' , true)) + +diag(A::LowerTriangular) = diag(A.data) +diag(A::UnitLowerTriangular) = ones(eltype(A), size(A,1)) +diag(A::UpperTriangular) = diag(A.data) +diag(A::UnitUpperTriangular) = ones(eltype(A), size(A,1)) # Unary operations --{T, S, UpLo}(A::Triangular{T, S, UpLo, false}) = Triangular{T, S, UpLo, false}(-A.data) -function -{T, S, UpLo}(A::Triangular{T, S, UpLo, true}) +-(A::LowerTriangular) = LowerTriangular(-A.data) +-(A::UpperTriangular) = UpperTriangular(-A.data) +function -(A::UnitLowerTriangular) Anew = -A.data for i = 1:size(A, 1) Anew[i, i] = -1 end - Triangular{T, S, UpLo, false}(Anew) + LowerTriangular(Anew) +end +function -(A::UnitUpperTriangular) + Anew = -A.data + for i = 1:size(A, 1) + Anew[i, i] = -1 + end + UpperTriangular(Anew) end # Binary operations -+{TA, TB, SA, SB, uplo}(A::Triangular{TA, SA, uplo, false}, B::Triangular{TB, SB, uplo, false}) = Triangular(A.data + B.data, uplo) -+{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, false}, B::Triangular{TB, SB, :U, true}) = Triangular(A.data + triu(B.data, 1) + I, :U) -+{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, false}, B::Triangular{TB, SB, :L, true}) = Triangular(A.data + tril(B.data, -1) + I, :L) -+{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, false}) = Triangular(triu(A.data, 1) + B.data + I, :U) -+{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, false}) = Triangular(tril(A.data, -1) + B.data + I, :L) -+{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, true}) = Triangular(triu(A.data, 1) + triu(B.data, 1) + 2I, :U) -+{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, true}) = Triangular(tril(A.data, -1) + tril(B.data, -1) + 2I, :L) -+{TA, SA, TB, SB, uplo1, uplo2, IsUnit1, IsUnit2}(A::Triangular{TA, SA, uplo1, IsUnit1}, B::Triangular{TB, SB, uplo2, IsUnit2}) = full(A) + full(B) --{TA, SA, TB, SB, uplo}(A::Triangular{TA, SA, uplo, false}, B::Triangular{TB, SB, uplo, false}) = Triangular(A.data - B.data, uplo) --{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, false}, B::Triangular{TB, SB, :U, true}) = Triangular(A.data - triu(B.data, 1) - I, :U) --{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, false}, B::Triangular{TB, SB, :L, true}) = Triangular(A.data - tril(B.data, -1) - I, :L) --{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, false}) = Triangular(triu(A.data, 1) - B.data + I, :U) --{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, false}) = Triangular(tril(A.data, -1) - B.data + I, :L) --{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, true}) = Triangular(triu(A.data, 1) - triu(B.data, 1), :U) --{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, true}) = Triangular(tril(A.data, -1) - tril(B.data, -1), :L) --{TA, SA, TB, SB, uplo1, uplo2, IsUnit1, IsUnit2}(A::Triangular{TA, SA, uplo1, IsUnit1}, B::Triangular{TB, SB, uplo2, IsUnit2}) = full(A) - full(B) ++(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data + B.data) ++(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(A.data + B.data) ++(A::UpperTriangular, B::UnitUpperTriangular) = UpperTriangular(A.data + triu(B.data, 1) + I) ++(A::LowerTriangular, B::UnitLowerTriangular) = LowerTriangular(A.data + tril(B.data, -1) + I) ++(A::UnitUpperTriangular, B::UpperTriangular) = UpperTriangular(triu(A.data, 1) + B.data + I) ++(A::UnitLowerTriangular, B::LowerTriangular) = LowerTriangular(tril(A.data, -1) + B.data + I) ++(A::UnitUpperTriangular, B::UnitUpperTriangular) = UpperTriangular(triu(A.data, 1) + triu(B.data, 1) + 2I) ++(A::UnitLowerTriangular, B::UnitLowerTriangular) = LowerTriangular(tril(A.data, -1) + tril(B.data, -1) + 2I) + +-(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data - B.data) +-(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(A.data - B.data) +-(A::UpperTriangular, B::UnitUpperTriangular) = UpperTriangular(A.data - triu(B.data, 1) - I) +-(A::LowerTriangular, B::UnitLowerTriangular) = LowerTriangular(A.data - tril(B.data, -1) - I) +-(A::UnitUpperTriangular, B::UpperTriangular) = UpperTriangular(triu(A.data, 1) - B.data + I) +-(A::UnitLowerTriangular, B::LowerTriangular) = LowerTriangular(tril(A.data, -1) - B.data + I) +-(A::UnitUpperTriangular, B::UnitUpperTriangular) = UpperTriangular(triu(A.data, 1) - triu(B.data, 1)) +-(A::UnitLowerTriangular, B::UnitLowerTriangular) = LowerTriangular(tril(A.data, -1) - tril(B.data, -1)) ###################### # BlasFloat routines # ###################### -for (uplos, uploc) in ((:(:U), 'U'), (:(:L), 'L')) - for (isunitb, isunitc) in ((true, 'U'), (false, 'N')) - @eval begin +A_mul_B!(A::Tridiagonal, B::TriangularUnion) = A*full!(B) +A_mul_B!(C::AbstractVecOrMat, A::TriangularUnion, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B)) +A_mul_Bc!(C::AbstractVecOrMat, A::TriangularUnion, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B)) + +for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), + (:UnitLowerTriangular, 'L', 'U'), + (:UpperTriangular, 'U', 'N'), + (:UnitUpperTriangular, 'U', 'U')) + @eval begin # Vector multiplication - A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, b::StridedVector{T}) = BLAS.trmv!($uploc, 'N', $isunitc, A.data, b) + A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, b::StridedVector{T}) = BLAS.trmv!($uploc, 'N', $isunitc, A.data, b) # Matrix multiplication - A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B) - A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A) + A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B) + A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A) - Ac_mul_B!{T<:BlasComplex,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B) - Ac_mul_B!{T<:BlasReal,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B) + Ac_mul_B!{T<:BlasComplex,S<:StridedMatrix}(A::$t{T,S}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B) + Ac_mul_B!{T<:BlasReal,S<:StridedMatrix}(A::$t{T,S}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B) - A_mul_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A) - A_mul_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A) + A_mul_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A) + A_mul_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A) # Left division - A_ldiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'N', $isunitc, A.data, B) - Ac_ldiv_B!{T<:BlasReal,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'T', $isunitc, A.data, B) - Ac_ldiv_B!{T<:BlasComplex,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'C', $isunitc, A.data, B) + A_ldiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'N', $isunitc, A.data, B) + Ac_ldiv_B!{T<:BlasReal,S<:StridedMatrix}(A::$t{T,S}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'T', $isunitc, A.data, B) + Ac_ldiv_B!{T<:BlasComplex,S<:StridedMatrix}(A::$t{T,S}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'C', $isunitc, A.data, B) # Right division - A_rdiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedVecOrMat{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trsm!('R', $uploc, 'N', $isunitc, one(T), B.data, A) - A_rdiv_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trsm!('R', $uploc, 'T', $isunitc, one(T), B.data, A) - A_rdiv_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trsm!('R', $uploc, 'C', $isunitc, one(T), B.data, A) + A_rdiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedVecOrMat{T}, B::$t{T,S}) = BLAS.trsm!('R', $uploc, 'N', $isunitc, one(T), B.data, A) + A_rdiv_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trsm!('R', $uploc, 'T', $isunitc, one(T), B.data, A) + A_rdiv_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trsm!('R', $uploc, 'C', $isunitc, one(T), B.data, A) # Matrix inverse - inv{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}) = Triangular{T,S,$uplos,$isunitb}(LAPACK.trtri!($uploc, $isunitc, copy(A.data))) + inv{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}) = $t{T,S}(LAPACK.trtri!($uploc, $isunitc, copy(A.data))) # Error bounds for triangular solve - errorbounds{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = LAPACK.trrfs!($uploc, 'N', $isunitc, A.data, B, X) + errorbounds{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = LAPACK.trrfs!($uploc, 'N', $isunitc, A.data, B, X) # Condition numbers - function cond{T<:BlasFloat,S}(A::Triangular{T,S,$uplos,$isunitb}, p::Real=2) + function cond{T<:BlasFloat,S}(A::$t{T,S}, p::Real=2) chksquare(A) if p==1 return inv(LAPACK.trcon!('O', $uploc, $isunitc, A.data)) @@ -146,68 +231,115 @@ for (uplos, uploc) in ((:(:U), 'U'), (:(:L), 'L')) end end end - end end -errorbounds{T<:Union(BigFloat, Complex{BigFloat}),S<:StridedMatrix}(A::Triangular{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = error("not implemented yet! Please submit a pull request.") -function errorbounds{TA<:Number,S<:StridedMatrix,UpLo,IsUnit,TX<:Number,TB<:Number}(A::Triangular{TA,S,UpLo,IsUnit}, X::StridedVecOrMat{TX}, B::StridedVecOrMat{TB}) +errorbounds{T<:Union(BigFloat, Complex{BigFloat}),S<:StridedMatrix}(A::TriangularUnion{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = error("not implemented yet! Please submit a pull request.") +function errorbounds{TA<:Number,S<:StridedMatrix,TX<:Number,TB<:Number}(A::TriangularUnion{TA,S}, X::StridedVecOrMat{TX}, B::StridedVecOrMat{TB}) TAXB = promote_type(TA, TB, TX, Float32) errorbounds(convert(AbstractMatrix{TAXB}, A), convert(AbstractArray{TAXB}, X), convert(AbstractArray{TAXB}, B)) end # Eigensystems ## Notice that trecv works for quasi-triangular matrices and therefore the lower sub diagonal must be zeroed before calling the subroutine -eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:U,false}) = LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data)) -eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:U,true}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data))) -eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:L,false}) = LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)') -eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:L,true}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)')) +eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::UpperTriangular{T,S}) = LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data)) +eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::UnitUpperTriangular{T,S}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data))) +eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::LowerTriangular{T,S}) = LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)') +eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::UnitLowerTriangular{T,S}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)')) #################### # Generic routines # #################### -(*){T,S,UpLo}(A::Triangular{T,S,UpLo,false}, x::Number) = Triangular(A.data*x, UpLo) -function (*){T,S,UpLo}(A::Triangular{T,S,UpLo,true}, x::Number) +(*)(A::UpperTriangular, x::Number) = UpperTriangular(A.data*x) +(*)(A::LowerTriangular, x::Number) = LowerTriangular(A.data*x) +function (*)(A::UnitUpperTriangular, x::Number) + B = A.data*x + for i = 1:size(A, 1) + B[i,i] = x + end + UpperTriangular(B) +end +function (*)(A::UnitLowerTriangular, x::Number) B = A.data*x for i = 1:size(A, 1) B[i,i] = x end - Triangular(B, UpLo) + LowerTriangular(B) +end +(*)(x::Number, A::UpperTriangular) = UpperTriangular(x*A.data) +(*)(x::Number, A::LowerTriangular) = LowerTriangular(x*A.data) +function (*)(x::Number, A::UnitUpperTriangular) + B = x*A.data + for i = 1:size(A, 1) + B[i,i] = x + end + UpperTriangular(B) end -(*){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,false}) = Triangular(x*A.data, UpLo) -function (*){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,true}) +function (*)(x::Number, A::UnitLowerTriangular) B = x*A.data for i = 1:size(A, 1) B[i,i] = x end - Triangular(B, UpLo) + LowerTriangular(B) end -(/){T,S,UpLo}(A::Triangular{T,S,UpLo,false}, x::Number) = Triangular(A.data/x, UpLo) -function (/){T,S,UpLo}(A::Triangular{T,S,UpLo,true}, x::Number) +(/)(A::UpperTriangular, x::Number) = UpperTriangular(A.data/x) +(/)(A::LowerTriangular, x::Number) = LowerTriangular(A.data/x) +function (/)(A::UnitUpperTriangular, x::Number) B = A.data*x invx = inv(x) for i = 1:size(A, 1) B[i,i] = x end - Triangular(B, UpLo) + UpperTriangular(B) end -(\){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,false}) = Triangular(x\A.data, UpLo) -function (\){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,true}) +function (/)(A::UnitLowerTriangular, x::Number) + B = A.data*x + invx = inv(x) + for i = 1:size(A, 1) + B[i,i] = x + end + LowerTriangular(B) +end +(\)(x::Number, A::UpperTriangular) = UpperTriangular(x\A.data) +(\)(x::Number, A::LowerTriangular) = LowerTriangular(x\A.data) +function (\)(x::Number, A::UnitUpperTriangular) B = x\A.data invx = inv(x) for i = 1:size(A, 1) B[i,i] = invx end - Triangular(B, UpLo) + UpperTriangular(B) +end +function (\)(x::Number, A::UnitLowerTriangular) + B = x\A.data + invx = inv(x) + for i = 1:size(A, 1) + B[i,i] = invx + end + LowerTriangular(B) end ## Generic triangular multiplication -function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::StridedVecOrMat) +function A_mul_B!(A::UpperTriangular, B::StridedVecOrMat) + m, n = size(B, 1), size(B, 2) + m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) + for j = 1:n + for i = 1:m + Bij = A.data[i,i]*B[i,j] + for k = i + 1:m + Bij += A.data[i,k]*B[k,j] + end + B[i,j] = Bij + end + end + B +end +function A_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) for j = 1:n for i = 1:m - Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j]) + Bij = B[i,j] for k = i + 1:m Bij += A.data[i,k]*B[k,j] end @@ -217,12 +349,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::S B end -function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::StridedVecOrMat) +function A_mul_B!(A::LowerTriangular, B::StridedVecOrMat) + m, n = size(B, 1), size(B, 2) + m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) + for j = 1:n + for i = m:-1:1 + Bij = A.data[i,i]*B[i,j] + for k = 1:i - 1 + Bij += A.data[i,k]*B[k,j] + end + B[i,j] = Bij + end + end + B +end +function A_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) for j = 1:n for i = m:-1:1 - Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j]) + Bij = B[i,j] for k = 1:i - 1 Bij += A.data[i,k]*B[k,j] end @@ -232,12 +378,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::S B end -function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::StridedVecOrMat) +function Ac_mul_B!(A::UpperTriangular, B::StridedVecOrMat) + m, n = size(B, 1), size(B, 2) + m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) + for j = 1:n + for i = m:-1:1 + Bij = A.data[i,i]*B[i,j] + for k = 1:i - 1 + Bij += A.data[k,i]'B[k,j] + end + B[i,j] = Bij + end + end + B +end +function Ac_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) for j = 1:n for i = m:-1:1 - Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j]) + Bij = B[i,j] for k = 1:i - 1 Bij += A.data[k,i]'B[k,j] end @@ -247,12 +407,26 @@ function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B:: B end -function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::StridedVecOrMat) +function Ac_mul_B!(A::LowerTriangular, B::StridedVecOrMat) + m, n = size(B, 1), size(B, 2) + m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) + for j = 1:n + for i = 1:m + Bij = A.data[i,i]*B[i,j] + for k = i + 1:m + Bij += A.data[k,i]'B[k,j] + end + B[i,j] = Bij + end + end + B +end +function Ac_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit")) for j = 1:n for i = 1:m - Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j]) + Bij = B[i,j] for k = i + 1:m Bij += A.data[k,i]'B[k,j] end @@ -262,12 +436,26 @@ function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B:: B end -function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit}) +function A_mul_B!(A::StridedMatrix, B::UpperTriangular) + m, n = size(A) + size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = n:-1:1 + Aij = A[i,j]*B[j,j] + for k = 1:j - 1 + Aij += A[i,k]*B.data[k,j] + end + A[i,j] = Aij + end + end + A +end +function A_mul_B!(A::StridedMatrix, B::UnitUpperTriangular) m, n = size(A) size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m for j = n:-1:1 - Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j]) + Aij = A[i,j] for k = 1:j - 1 Aij += A[i,k]*B.data[k,j] end @@ -277,12 +465,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T, A end -function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit}) +function A_mul_B!(A::StridedMatrix, B::LowerTriangular) + m, n = size(A) + size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = 1:n + Aij = A[i,j]*B[j,j] + for k = j + 1:n + Aij += A[i,k]*B.data[k,j] + end + A[i,j] = Aij + end + end + A +end +function A_mul_B!(A::StridedMatrix, B::UnitLowerTriangular) m, n = size(A) size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m for j = 1:n - Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j]) + Aij = A[i,j] for k = j + 1:n Aij += A[i,k]*B.data[k,j] end @@ -292,12 +494,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T, A end -function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit}) +function A_mul_Bc!(A::StridedMatrix, B::UpperTriangular) m, n = size(A) size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m for j = 1:n - Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j]) + Aij = A[i,j]*B[j,j] + for k = j + 1:n + Aij += A[i,k]*B.data[j,k]' + end + A[i,j] = Aij + end + end + A +end +function A_mul_Bc!(A::StridedMatrix, B::UnitUpperTriangular) + m, n = size(A) + size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = 1:n + Aij = A[i,j] for k = j + 1:n Aij += A[i,k]*B.data[j,k]' end @@ -307,12 +523,26 @@ function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T A end -function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit}) +function A_mul_Bc!(A::StridedMatrix, B::LowerTriangular) m, n = size(A) size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m for j = n:-1:1 - Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j]) + Aij = A[i,j]*B[j,j] + for k = 1:j - 1 + Aij += A[i,k]*B.data[j,k]' + end + A[i,j] = Aij + end + end + A +end +function A_mul_Bc!(A::StridedMatrix, B::UnitLowerTriangular) + m, n = size(A) + size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = n:-1:1 + Aij = A[i,j] for k = 1:j - 1 Aij += A[i,k]*B.data[j,k]' end @@ -321,43 +551,70 @@ function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T end A end - -A_mul_B!(A::Tridiagonal, B::Triangular) = A*full!(B) -A_mul_B!(C::AbstractVecOrMat, A::Triangular, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B)) -A_mul_Bc!(C::AbstractVecOrMat, A::Triangular, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B)) #Generic solver using naive substitution -function naivesub!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::AbstractVector, x::AbstractVector=b) +function naivesub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector=b) N = size(A, 2) - N==length(b)==length(x) || throw(DimensionMismatch()) - - if UpLo == :L #do forward substitution - for j = 1:N - x[j] = b[j] - for k = 1:j-1 - x[j] -= A[j,k] * x[k] - end - if !IsUnit - x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j] - end + N == length(b) == length(x) || throw(DimensionMismatch()) + for j = N:-1:1 + x[j] = b[j] + for k = j+1:1:N + x[j] -= A[j,k] * x[k] end - elseif UpLo == :U #do backward substitution - for j = N:-1:1 - x[j] = b[j] - for k = j+1:1:N - x[j] -= A[j,k] * x[k] - end - if !IsUnit - x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j] - end + x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j] + end + x +end +function naivesub!(A::UnitUpperTriangular, b::AbstractVector, x::AbstractVector=b) + N = size(A, 2) + N == length(b) == length(x) || throw(DimensionMismatch()) + for j = N:-1:1 + x[j] = b[j] + for k = j+1:1:N + x[j] -= A[j,k] * x[k] + end + end + x +end +function naivesub!(A::LowerTriangular, b::AbstractVector, x::AbstractVector=b) + N = size(A, 2) + N == length(b) == length(x) || throw(DimensionMismatch()) + for j = 1:N + x[j] = b[j] + for k = 1:j-1 + x[j] -= A[j,k] * x[k] + end + x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j] + end + x +end +function naivesub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector=b) + N = size(A, 2) + N == length(b) == length(x) || throw(DimensionMismatch()) + for j = 1:N + x[j] = b[j] + for k = 1:j-1 + x[j] -= A[j,k] * x[k] end - else - throw(ArgumentError("Unknown UpLo=$(UpLo)")) end x end -function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit}) +function A_rdiv_B!(A::StridedMatrix, B::UpperTriangular) + m, n = size(A) + size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = 1:n + Aij = A[i,j] + for k = 1:j - 1 + Aij -= A[i,k]*B.data[k,j] + end + A[i,j] = Aij/B[j,j] + end + end + A +end +function A_rdiv_B!(A::StridedMatrix, B::UnitUpperTriangular) m, n = size(A) size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m @@ -366,13 +623,13 @@ function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit}) for k = 1:j - 1 Aij -= A[i,k]*B.data[k,j] end - A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j]) + A[i,j] = Aij end end A end -function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit}) +function A_rdiv_B!(A::StridedMatrix, B::LowerTriangular) m, n = size(A) size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m @@ -381,13 +638,41 @@ function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit}) for k = j + 1:n Aij -= A[i,k]*B.data[k,j] end - A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j]) + A[i,j] = Aij/B[j,j] + end + end + A +end +function A_rdiv_B!(A::StridedMatrix, B::UnitLowerTriangular) + m, n = size(A) + size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = n:-1:1 + Aij = A[i,j] + for k = j + 1:n + Aij -= A[i,k]*B.data[k,j] + end + A[i,j] = Aij end end A end -function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit}) +function A_rdiv_Bc!(A::StridedMatrix, B::UpperTriangular) + m, n = size(A) + size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = n:-1:1 + Aij = A[i,j] + for k = j + 1:n + Aij -= A[i,k]*B.data[j,k]' + end + A[i,j] = Aij/B[j,j] + end + end + A +end +function A_rdiv_Bc!(A::StridedMatrix, B::UnitUpperTriangular) m, n = size(A) size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m @@ -396,13 +681,27 @@ function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit}) for k = j + 1:n Aij -= A[i,k]*B.data[j,k]' end - A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j]) + A[i,j] = Aij end end A end -function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit}) +function A_rdiv_Bc!(A::StridedMatrix, B::LowerTriangular) + m, n = size(A) + size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) + for i = 1:m + for j = 1:n + Aij = A[i,j] + for k = 1:j - 1 + Aij -= A[i,k]*B.data[j,k]' + end + A[i,j] = Aij/B[j,j] + end + end + A +end +function A_rdiv_Bc!(A::StridedMatrix, B::UnitLowerTriangular) m, n = size(A) size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit")) for i = 1:m @@ -411,25 +710,30 @@ function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit}) for k = 1:j - 1 Aij -= A[i,k]*B.data[j,k]' end - A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j]) + A[i,j] = Aij end end A end # Promotion -## Promotion methods in matmul doesn't apply to triangular multiplication since it is inplace. Hence we have to make very similar definitions, but without allocation of a result array. For multiplication and unit diagonal division the element type doesn't have to be stable under division whereas that is necessary in the general triangular solve problem. +## Promotion methods in matmul don't apply to triangular multiplication since it is inplace. Hence we have to make very similar definitions, but without allocation of a result array. For multiplication and unit diagonal division the element type doesn't have to be stable under division whereas that is necessary in the general triangular solve problem. ## Some Triangular-Triangular cases. We might want to write taylored methods for these cases, but I'm not sure it is worth it. -*(A::Tridiagonal, B::Triangular) = A_mul_B!(full(A), B) +for t in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular) + @eval begin + *(A::Tridiagonal, B::$t) = A_mul_B!(full(A), B) + end +end + for f in (:*, :Ac_mul_B, :At_mul_B, :A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :\, :Ac_ldiv_B, :At_ldiv_B) @eval begin - ($f)(A::Triangular, B::Triangular) = ($f)(A, full(B)) + ($f)(A::TriangularUnion, B::TriangularUnion) = ($f)(A, full(B)) end end for f in (:A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :/, :A_rdiv_Bc, :A_rdiv_Bt) @eval begin - ($f)(A::Triangular, B::Triangular) = ($f)(full(A), B) + ($f)(A::TriangularUnion, B::TriangularUnion) = ($f)(full(A), B) end end @@ -437,7 +741,7 @@ end ### Multiplication with triangle to the left and hence rhs cannot be transposed. for (f, g) in ((:*, :A_mul_B!), (:Ac_mul_B, :Ac_mul_B!), (:At_mul_B, :At_mul_B!)) @eval begin - function ($f){TA,TB}(A::Triangular{TA}, B::StridedVecOrMat{TB}) + function ($f){TA,TB}(A::TriangularUnion{TA}, B::StridedVecOrMat{TB}) TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -446,7 +750,15 @@ end ### Left division with triangle to the left hence rhs cannot be transposed. No quotients. for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!)) @eval begin - function ($f){TA,TB,S,UpLo}(A::Triangular{TA,S,UpLo,true}, B::StridedVecOrMat{TB}) + function ($f){TA,TB,S}(A::UnitUpperTriangular{TA,S}, B::StridedVecOrMat{TB}) + TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) + ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) + end + end +end +for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!)) + @eval begin + function ($f){TA,TB,S}(A::UnitLowerTriangular{TA,S}, B::StridedVecOrMat{TB}) TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -455,7 +767,15 @@ end ### Left division with triangle to the left hence rhs cannot be transposed. Quotients. for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!)) @eval begin - function ($f){TA,TB,S,UpLo}(A::Triangular{TA,S,UpLo,false}, B::StridedVecOrMat{TB}) + function ($f){TA,TB,S}(A::UpperTriangular{TA,S}, B::StridedVecOrMat{TB}) + TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA)) + ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) + end + end +end +for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!)) + @eval begin + function ($f){TA,TB,S}(A::LowerTriangular{TA,S}, B::StridedVecOrMat{TB}) TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -464,7 +784,7 @@ end ### Multiplication with triangle to the rigth and hence lhs cannot be transposed. for (f, g) in ((:*, :A_mul_B!), (:A_mul_Bc, :A_mul_Bc!), (:A_mul_Bt, :A_mul_Bt!)) @eval begin - function ($f){TA,TB}(A::StridedVecOrMat{TA}, B::Triangular{TB}) + function ($f){TA,TB}(A::StridedVecOrMat{TA}, B::TriangularUnion{TB}) TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -473,7 +793,15 @@ end ### Right division with triangle to the right hence lhs cannot be transposed. No quotients. for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!)) @eval begin - function ($f){TA,TB,S,UpLo}(A::StridedVecOrMat{TA}, B::Triangular{TB,S,UpLo,true}) + function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::UnitUpperTriangular{TB,S}) + TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) + ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) + end + end +end +for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!)) + @eval begin + function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::UnitLowerTriangular{TB,S}) TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -482,52 +810,76 @@ end ### Right division with triangle to the right hence lhs cannot be transposed. Quotients. for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!)) @eval begin - function ($f){TA,TB,S,UpLo}(A::StridedVecOrMat{TA}, B::Triangular{TB,S,UpLo,false}) + function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::UpperTriangular{TB,S}) + TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA)) + ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) + end + end +end +for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!)) + @eval begin + function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::LowerTriangular{TB,S}) TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end end end -function sqrtm{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) +function sqrtm{T}(A::UpperTriangular{T}) n = size(A, 1) TT = typeof(sqrt(zero(T))) R = zeros(TT, n, n) - if UpLo == :U - for j = 1:n - (T<:Complex || A[j,j]>=0) ? (R[j,j] = IsUnit ? one(T) : sqrt(A[j,j])) : throw(SingularException(j)) - for i = j-1:-1:1 - r = A[i,j] - for k = i+1:j-1 - r -= R[i,k]*R[k,j] - end - r==0 || (R[i,j] = r / (R[i,i] + R[j,j])) + for j = 1:n + R[j,j] = sqrt(A[j,j]) + for i = j-1:-1:1 + r = A[i,j] + for k = i+1:j-1 + r -= R[i,k]*R[k,j] + end + r==0 || (R[i,j] = r / (R[i,i] + R[j,j])) + end + end + return UpperTriangular(R) +end +function sqrtm{T}(A::UnitUpperTriangular{T}) + n = size(A, 1) + TT = typeof(sqrt(zero(T))) + R = zeros(TT, n, n) + for j = 1:n + R[j,j] = one(T) + for i = j-1:-1:1 + r = A[i,j] + for k = i+1:j-1 + r -= R[i,k]*R[k,j] end + r==0 || (R[i,j] = r / (R[i,i] + R[j,j])) end - return Triangular(R, :U, IsUnit) - else # UpLo == :L #Not the usual case - return sqrtm(A.').' end + return UnitUpperTriangular(R) end +sqrtm(A::LowerTriangular) = sqrtm(A.').' +sqrtm(A::UnitLowerTriangular) = sqrtm(A.').' #Generic eigensystems -eigvals(A::Triangular) = diag(A) -function eigvecs{T}(A::Triangular{T}) +eigvals(A::TriangularUnion) = diag(A) +function eigvecs{T}(A::TriangularUnion{T}) TT = promote_type(T, Float32) TT <: BlasFloat && return eigvecs(convert(AbstractMatrix{TT}, A)) error("type not handled yet. Please submit a pull request.") end -det{T,S,UpLo}(A::Triangular{T,S,UpLo,true}) = one(T)*one(T) -det{T,S,UpLo}(A::Triangular{T,S,UpLo,false}) = prod(diag(A.data)) +det{T}(A::UnitUpperTriangular{T}) = one(T)*one(T) +det{T}(A::UnitLowerTriangular{T}) = one(T)*one(T) +det{T}(A::UpperTriangular{T}) = prod(diag(A.data)) +det{T}(A::LowerTriangular{T}) = prod(diag(A.data)) -eigfact(A::Triangular) = Eigen(eigvals(A), eigvecs(A)) +eigfact(A::TriangularUnion) = Eigen(eigvals(A), eigvecs(A)) #Generic singular systems for func in (:svd, :svdfact, :svdfact!, :svdvals) @eval begin - ($func)(A::Triangular) = ($func)(full(A)) + ($func)(A::TriangularUnion) = ($func)(full(A)) end end -factorize(A::Triangular) = A +factorize(A::TriangularUnion) = A diff --git a/base/linalg/uniformscaling.jl b/base/linalg/uniformscaling.jl index a13a2f0d42d51..4bf9627540511 100644 --- a/base/linalg/uniformscaling.jl +++ b/base/linalg/uniformscaling.jl @@ -84,7 +84,7 @@ inv(J::UniformScaling) = UniformScaling(inv(J.λ)) /(J::UniformScaling, x::Number) = UniformScaling(J.λ/x) \(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ) -\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), J::UniformScaling) = inv(A)*J.λ +\{T<:Number}(A::Union(Bidiagonal{T},TriangularUnion{T}), J::UniformScaling) = inv(A)*J.λ \(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A \(A::AbstractMatrix, J::UniformScaling) = inv(A)*J.λ diff --git a/doc/manual/linear-algebra.rst b/doc/manual/linear-algebra.rst index 9e6dd72f4bec0..733dc96eca888 100644 --- a/doc/manual/linear-algebra.rst +++ b/doc/manual/linear-algebra.rst @@ -46,7 +46,8 @@ for them in LAPACK are available. ======================= ================================================================================== :class:`Hermitian` `Hermitian matrix `_ -:class:`Triangular` Upper/lower `triangular matrix `_ +:class:`UpperTriangular Upper `triangular matrix `_ +:class:`LowerTriangular Lower `triangular matrix `_ :class:`Tridiagonal` `Tridiagonal matrix `_ :class:`SymTridiagonal` Symmetric tridiagonal matrix :class:`Bidiagonal` Upper/lower `bidiagonal matrix `_ @@ -64,7 +65,9 @@ Elementary operations | :class:`Hermitian` | | | | MV | :func:`inv`, | | | | | | | :func:`sqrtm`, :func:`expm` | +-------------------------+-------+-------+-------+-------+--------------------------------+ -| :class:`Triangular` | | | MV | MV | :func:`inv`, :func:`det` | +| :class:`UpperTriangular`| | | MV | MV | :func:`inv`, :func:`det` | ++-------------------------+-------+-------+-------+-------+--------------------------------+ +| :class:`LowerTriangular`| | | MV | MV | :func:`inv`, :func:`det` | +-------------------------+-------+-------+-------+-------+--------------------------------+ | :class:`SymTridiagonal` | M | M | MS | MV | :func:`eigmax`, :func:`eigmin` | +-------------------------+-------+-------+-------+-------+--------------------------------+ @@ -96,7 +99,9 @@ Matrix factorizations +=========================+========+=============+=================+=================+=============+=================+ | :class:`Hermitian` | HE | | ARI | | | | +-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+ -| :class:`Triangular` | TR | | | | | | +| :class:`UpperTriangular`| TR | A | A | A | | | ++-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+ +| :class:`LowerTriangular`| TR | A | A | A | | | +-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+ | :class:`SymTridiagonal` | ST | A | ARI | AV | | | +-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+ diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 4f1b83831787a..e118da5312267 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -39,7 +39,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f .. function:: factorize(A) - Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, Triangular) of A, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example: ``A=factorize(A); x=A\\b; y=A\\C``. + Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, LowerTriangular, UpperTriangular) of A, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example: ``A=factorize(A); x=A\\b; y=A\\C``. .. function:: factorize!(A) diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl index bbb8be43c1734..785a6c3dbb16a 100644 --- a/test/linalg/triangular.jl +++ b/test/linalg/triangular.jl @@ -1,6 +1,6 @@ debug = false using Base.Test -using Base.LinAlg: BlasFloat, errorbounds, full!, naivesub!, transpose! +using Base.LinAlg: BlasFloat, errorbounds, full!, naivesub!, transpose!, UnitUpperTriangular, UnitLowerTriangular debug && println("Triangular matrices") @@ -8,145 +8,147 @@ n = 9 srand(123) debug && println("Test basic type functionality") -@test_throws DimensionMismatch Triangular(randn(5, 4), :L) -@test Triangular(randn(3, 3), :L) |> t -> [size(t, i) for i = 1:3] == [size(full(t), i) for i = 1:3] +@test_throws DimensionMismatch LowerTriangular(randn(5, 4)) +@test LowerTriangular(randn(3, 3)) |> t -> [size(t, i) for i = 1:3] == [size(full(t), i) for i = 1:3] # The following test block tries to call all methods in base/linalg/triangular.jl in order for a combination of input element types. Keep the ordering when adding code. for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) for elty2 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) - for uplo1 in (:U, :L) - for uplo2 in (:U, :L) - for isunit1 in (true, false) - for isunit2 in (true, false) - A1 = Triangular(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, uplo1)), uplo1, isunit1) - A2 = Triangular(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, uplo2)), uplo2, isunit2) - - for eltyB in (Float32, Float64, Complex64, Complex128) - B = convert(Matrix{eltyB}, elty1 <: Complex ? real(A1)*ones(n, n) : A1*ones(n, n)) - - debug && println("A1: $elty1, $uplo1, $isunit1.A2: $elty2, $uplo2, $isunit2. B: $eltyB") - - # Convert - @test convert(Triangular{elty1}, A1) == A1 - if elty1 <: Real && !(elty2 <: Integer) - @test convert(Triangular{elty2}, A1) == Triangular(convert(Matrix{elty2}, A1.data), uplo1, isunit1) - elseif elty2 <: Real && !(elty1 <: Integer) - @test_throws InexactError convert(Triangular{elty2}, A1) == Triangular(convert(Matrix{elty2}, A1.data), uplo1, isunit1) - end - @test convert(Matrix, A1) == full(A1) - - # full! - @test full!(copy(A1)) == full(A1) - - # fill! - @test full!(fill!(copy(A1), 1)) == full(Triangular(ones(size(A1)...), uplo1, isunit1)) - - # similar - @test isa(similar(A1), Triangular{elty1, Matrix{elty1}, uplo1, isunit1}) - - # Linear indexing - for i = 1:length(A1) - @test A1[i] == full(A1)[i] - end - - # istril/istriu - if uplo1 == :L - @test istril(A1) - @test !istriu(A1) - else - @test istriu(A1) - @test !istril(A1) - end - - # (c)transpose - @test full(A1') == full(A1)' - @test full(A1.') == full(A1).' - @test transpose!(copy(A1)) == A1.' - - # diag - @test diag(A1) == diag(full(A1)) - - # real - @test full(real(A1)) == real(full(A1)) - - # Unary operations - @test full(-A1) == -full(A1) - - # Binary operations - @test full(A1 + A2) == full(A1) + full(A2) - @test full(A1 - A2) == full(A1) - full(A2) - @test A1*0.5 == full(A1*0.5) - @test 0.5*A1 == full(0.5*A1) - @test A1/0.5 == full(A1/0.5) - @test 0.5\A1 == full(0.5\A1) - - # Triangular-Triangular multiplication and division - elty1 != BigFloat && elty2 != BigFloat && @test_approx_eq full(A1*A2) full(A1)*full(A2) - @test_approx_eq full(A1'A2) full(A1)'full(A2) - @test_approx_eq full(A1*A2') full(A1)*full(A2)' - @test_approx_eq full(A1'A2') full(A1)'full(A2)' - @test_approx_eq full(A1/A2) full(A1)/full(A2) - - # Triangular-dense Matrix/vector multiplication - @test_approx_eq A1*B[:,1] full(A1)*B[:,1] - @test_approx_eq A1*B full(A1)*B - @test_approx_eq A1'B[:,1] full(A1)'B[:,1] - @test_approx_eq A1'B full(A1)'B - @test_approx_eq A1*B' full(A1)*B' - @test_approx_eq B*A1 B*full(A1) - @test_approx_eq B[:,1]'A1 B[:,1]'full(A1) - @test_approx_eq B'A1 B'full(A1) - @test_approx_eq B*A1' B*full(A1)' - @test_approx_eq B[:,1]'A1' B[:,1]'full(A1)' - @test_approx_eq B'A1' B'full(A1)' - - # ... and division - @test_approx_eq A1\B[:,1] full(A1)\B[:,1] - @test_approx_eq A1\B full(A1)\B - @test_approx_eq A1'\B[:,1] full(A1)'\B[:,1] - @test_approx_eq A1'\B full(A1)'\B - @test_approx_eq A1\B' full(A1)\B' - @test_approx_eq A1'\B' full(A1)'\B' - @test_approx_eq B/A1 B/full(A1) - @test_approx_eq B/A1' B/full(A1)' - @test_approx_eq B'/A1 B'/full(A1) - @test_approx_eq B'/A1' B'/full(A1)' - - # inversion - @test_approx_eq inv(A1) inv(lufact(full(A1))) - - # Error bounds - elty1 != BigFloat && errorbounds(A1, A1\B, B) - - # Determinant - @test_approx_eq det(A1) det(lufact(full(A1))) - - # Matri square root - @test_approx_eq sqrtm(A1) |> t->t*t A1 - - # eigenproblems - if elty1 != BigFloat # Not handled yet - vals, vecs = eig(A1) - if !isunit1 && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues. - @test_approx_eq_eps vecs*diagm(vals)/vecs full(A1) sqrt(eps(float(real(one(vals[1])))))*(norm(A1, Inf)*n)^2 - end - end - - # Condition number tests - can be VERY approximate - if elty1 <:BlasFloat - for p in (1.0, Inf) - @test_approx_eq_eps cond(A1, p) cond(A1, p) (cond(A1, p) + cond(A1, p)) - end - end - - if elty1 != BigFloat # Not implemented yet - svd(A1) - svdfact(A1) - elty1 <: BlasFloat && svdfact!(copy(A1)) - svdvals(A1) - end + for (t1, uplo1) in ((UpperTriangular, :U), + (UnitUpperTriangular, :U), + (LowerTriangular, :L), + (UnitLowerTriangular, :L)) + for (t2, uplo2) in ((UpperTriangular, :U), + (UnitUpperTriangular, :U), + (LowerTriangular, :L), + (UnitLowerTriangular, :L)) + A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't, uplo1))) + A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, uplo2))) + + for eltyB in (Float32, Float64, Complex64, Complex128) + B = convert(Matrix{eltyB}, elty1 <: Complex ? real(A1)*ones(n, n) : A1*ones(n, n)) + + debug && println("elty1: $elty1, A1: $t1, elty2: $elty2, A2: $t2, B: $eltyB") + + # Convert + @test convert(AbstractMatrix{elty1}, A1) == A1 + if elty1 <: Real && !(elty2 <: Integer) + @test convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data)) + elseif elty2 <: Real && !(elty1 <: Integer) + @test_throws InexactError convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data)) + end + @test convert(Matrix, A1) == full(A1) + + # full! + @test full!(copy(A1)) == full(A1) + + # fill! + @test full!(fill!(copy(A1), 1)) == full(t1(ones(size(A1)...))) + + # similar + @test isa(similar(A1), t1) + + # Linear indexing + for i = 1:length(A1) + @test A1[i] == full(A1)[i] + end + + # istril/istriu + if uplo1 == :L + @test istril(A1) + @test !istriu(A1) + else + @test istriu(A1) + @test !istril(A1) + end + + # (c)transpose + @test full(A1') == full(A1)' + @test full(A1.') == full(A1).' + @test transpose!(copy(A1)) == A1.' + + # diag + @test diag(A1) == diag(full(A1)) + + # real + @test full(real(A1)) == real(full(A1)) + + # Unary operations + @test full(-A1) == -full(A1) + + # Binary operations + @test full(A1 + A2) == full(A1) + full(A2) + @test full(A1 - A2) == full(A1) - full(A2) + @test A1*0.5 == full(A1*0.5) + @test 0.5*A1 == full(0.5*A1) + @test A1/0.5 == full(A1/0.5) + @test 0.5\A1 == full(0.5\A1) + + # Triangular-Triangular multiplication and division + elty1 != BigFloat && elty2 != BigFloat && @test_approx_eq full(A1*A2) full(A1)*full(A2) + @test_approx_eq full(A1'A2) full(A1)'full(A2) + @test_approx_eq full(A1*A2') full(A1)*full(A2)' + @test_approx_eq full(A1'A2') full(A1)'full(A2)' + @test_approx_eq full(A1/A2) full(A1)/full(A2) + + # Triangular-dense Matrix/vector multiplication + @test_approx_eq A1*B[:,1] full(A1)*B[:,1] + @test_approx_eq A1*B full(A1)*B + @test_approx_eq A1'B[:,1] full(A1)'B[:,1] + @test_approx_eq A1'B full(A1)'B + @test_approx_eq A1*B' full(A1)*B' + @test_approx_eq B*A1 B*full(A1) + @test_approx_eq B[:,1]'A1 B[:,1]'full(A1) + @test_approx_eq B'A1 B'full(A1) + @test_approx_eq B*A1' B*full(A1)' + @test_approx_eq B[:,1]'A1' B[:,1]'full(A1)' + @test_approx_eq B'A1' B'full(A1)' + + # ... and division + @test_approx_eq A1\B[:,1] full(A1)\B[:,1] + @test_approx_eq A1\B full(A1)\B + @test_approx_eq A1'\B[:,1] full(A1)'\B[:,1] + @test_approx_eq A1'\B full(A1)'\B + @test_approx_eq A1\B' full(A1)\B' + @test_approx_eq A1'\B' full(A1)'\B' + @test_approx_eq B/A1 B/full(A1) + @test_approx_eq B/A1' B/full(A1)' + @test_approx_eq B'/A1 B'/full(A1) + @test_approx_eq B'/A1' B'/full(A1)' + + # inversion + @test_approx_eq inv(A1) inv(lufact(full(A1))) + + # Error bounds + elty1 != BigFloat && errorbounds(A1, A1\B, B) + + # Determinant + @test_approx_eq_eps det(A1) det(lufact(full(A1))) sqrt(eps(real(float(one(elty1)))))*n*n + + # Matrix square root + @test_approx_eq sqrtm(A1) |> t->t*t A1 + + # eigenproblems + if elty1 != BigFloat # Not handled yet + vals, vecs = eig(A1) + if (t1 == UpperTriangular || t1 == LowerTriangular) && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues. + @test_approx_eq_eps vecs*diagm(vals)/vecs full(A1) sqrt(eps(float(real(one(vals[1])))))*(norm(A1, Inf)*n)^2 end end + + # Condition number tests - can be VERY approximate + if elty1 <:BlasFloat + for p in (1.0, Inf) + @test_approx_eq_eps cond(A1, p) cond(A1, p) (cond(A1, p) + cond(A1, p)) + end + end + + if elty1 != BigFloat # Not implemented yet + svd(A1) + svdfact(A1) + elty1 <: BlasFloat && svdfact!(copy(A1)) + svdvals(A1) + end end end end @@ -170,14 +172,14 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) debug && println("\ntype of A: ", eltya, " type of b: ", eltyb, "\n") debug && println("Solve upper triangular system") - Atri = Triangular(lufact(A)[:U], :U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned + Atri = UpperTriangular(lufact(A)[:U]) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned b = convert(Matrix{eltyb}, eltya <: Complex ? full(Atri)*ones(n, 2) : full(Atri)*ones(n, 2)) x = full(Atri) \ b debug && println("Test error estimates") if eltya != BigFloat && eltyb != BigFloat for i = 1:2 - @test norm(x[:,1] .- 1) <= errorbounds(Triangular(A, :U), x, b)[1][i] + @test norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i] end end debug && println("Test forward error [JIN 5705] if this is not a BigFloat") @@ -198,14 +200,14 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) end debug && println("Solve lower triangular system") - Atri = Triangular(lufact(A)[:U], :U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned + Atri = UpperTriangular(lufact(A)[:U]) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned b = convert(Matrix{eltyb}, eltya <: Complex ? full(Atri)*ones(n, 2) : full(Atri)*ones(n, 2)) x = full(Atri)\b debug && println("Test error estimates") if eltya != BigFloat && eltyb != BigFloat for i = 1:2 - @test norm(x[:,1] .- 1) <= errorbounds(Triangular(A, :U), x, b)[1][i] + @test norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i] end end diff --git a/test/linalg2.jl b/test/linalg2.jl index 3d11e4fd8bc9a..2e7a19bf2dd2a 100644 --- a/test/linalg2.jl +++ b/test/linalg2.jl @@ -403,7 +403,7 @@ S = sprandn(3,3,0.5) @test A/I == A @test I/λ === UniformScaling(1/λ) @test I\J === J -T = Triangular(randn(3,3),:L) +T = LowerTriangular(randn(3,3)) @test T\I == inv(T) @test I\A == A @test A\I == inv(A) diff --git a/test/linalg4.jl b/test/linalg4.jl index 29187946f5aa8..acab808152d18 100644 --- a/test/linalg4.jl +++ b/test/linalg4.jl @@ -154,7 +154,7 @@ end debug && println("Test interconversion between special matrix types") a=[1.0:n] A=Diagonal(a) -for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] +for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, LowerTriangular, UpperTriangular, Matrix] debug && println("newtype is $(newtype)") @test full(convert(newtype, A)) == full(A) end @@ -162,29 +162,29 @@ end for isupper in (true, false) debug && println("isupper is $(isupper)") A=Bidiagonal(a, [1.0:n-1], isupper) - for newtype in [Bidiagonal, Tridiagonal, Triangular, Matrix] + for newtype in [Bidiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix] debug && println("newtype is $(newtype)") @test full(convert(newtype, A)) == full(A) end A=Bidiagonal(a, zeros(n-1), isupper) #morally Diagonal - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix] debug && println("newtype is $(newtype)") @test full(convert(newtype, A)) == full(A) end end -A=SymTridiagonal(a, [1.0:n-1]) +A = SymTridiagonal(a, [1.0:n-1]) for newtype in [Tridiagonal, Matrix] @test full(convert(newtype, A)) == full(A) end -A=Tridiagonal(zeros(n-1), [1.0:n], zeros(n-1)) #morally Diagonal +A = Tridiagonal(zeros(n-1), [1.0:n], zeros(n-1)) #morally Diagonal for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix] @test full(convert(newtype, A)) == full(A) end -A=Triangular(full(Diagonal(a)), :L) #morally Diagonal -for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] +A = LowerTriangular(full(Diagonal(a))) #morally Diagonal +for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix] @test full(convert(newtype, A)) == full(A) end