Parametrize Triangular on matrix type

NEWS.md

@@ -260,6 +260,8 @@ Library improvements
       * new LAPACK wrappers
         - condition number estimate `cond(A::Triangular)` ([#5255])
 
+      * parametrize Triangular on matrix type ([#7064])
+
     * Dense linear algebra for generic matrix element types
 
       * LU factorization ([#5381] and [#5430])

base/array.jl

@@ -168,7 +168,7 @@ for (fname, felt) in ((:zeros,:zero), (:ones,:one))
     @eval begin
         ($fname){T}(::Type{T}, dims...)  = fill!(Array(T, dims...), ($felt)(T))
         ($fname)(dims...)                = fill!(Array(Float64, dims...), ($felt)(Float64))
-        ($fname){T}(x::AbstractArray{T}) = ($fname)(T, size(x))
+        ($fname){T}(A::AbstractArray{T}) = fill!(similar(A), ($felt)(T))
     end
 end
 

base/deprecated.jl

@@ -198,6 +198,8 @@ end
 @deprecate (/)(x::Number,A::Array)          x ./ A
 @deprecate (\)(A::Array,x::Number)          A .\ x
 
+@deprecate Triangular(A::Matrix)    Triangular(A, istril(A) ? :L : (istriu(A) ? istriu(A) : throw(ArgumentError("Matrix is not triangular"))))
+
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)
 deprecated_ls(args::String...) = run(`ls -l $args`)

base/linalg/bitarray.jl

@@ -40,7 +40,7 @@ end
 
 #aCb{T, S}(A::BitMatrix{T}, B::BitMatrix{S}) = aTb(A, B)
 
-function triu(B::BitMatrix, k::Integer)
+function triu(B::BitMatrix, k::Integer=0)
     m,n = size(B)
     A = falses(m,n)
     Ac = A.chunks
@@ -52,7 +52,7 @@ function triu(B::BitMatrix, k::Integer)
     A
 end
 
-function tril(B::BitMatrix, k::Integer)
+function tril(B::BitMatrix, k::Integer=0)
     m,n = size(B)
     A = falses(m, n)
     Ac = A.chunks

base/linalg/blas.jl

@@ -400,7 +400,7 @@ for (fname, elty) in ((:dtrmv_,:Float64),
                 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt},
                  Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                  &uplo, &trans, &diag, &n, 
-                 A, &max(1,stride(A,2)), x, &1)
+                 A, &max(1,stride(A,2)), x, &max(1,stride(x, 1)))
             x
         end
         function trmv(uplo::Char, trans::Char, diag::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty})

base/linalg/dense.jl

@@ -300,15 +300,15 @@ 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])))
+    R = full(sqrtm(Triangular(SchurF[:T], :U, false)))
     retmat = SchurF[:vectors]*R*SchurF[:vectors]'
     all(imag(retmat) .== 0) ? real(retmat) : retmat
 end
 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])))
+    R = full(sqrtm(Triangular(SchurF[:T], :U, false)))
     SchurF[:vectors]*R*SchurF[:vectors]'
 end
 sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)

base/linalg/diagonal.jl

@@ -8,6 +8,7 @@ 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)
 
 function similar{T}(D::Diagonal, ::Type{T}, d::(Int,Int))
     d[1] == d[2] || throw(ArgumentError("Diagonal matrix must be square"))
@@ -19,6 +20,8 @@ copy!(D1::Diagonal, D2::Diagonal) = (copy!(D1.diag, D2.diag); D1)
 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)
 
+fill!(D::Diagonal, x) = (fill!(D.diag, x); D)
+
 full(D::Diagonal) = diagm(D.diag)
 getindex(D::Diagonal, i::Integer, j::Integer) = i == j ? D.diag[i] : zero(eltype(D.diag))
 
@@ -33,6 +36,9 @@ ishermitian(D::Diagonal) = true
 issym(D::Diagonal) = true
 isposdef(D::Diagonal) = all(D.diag .> 0)
 
+tril!(D::Diagonal,i::Integer) = i == 0 ? D : zeros(D)
+triu!(D::Diagonal,i::Integer) = i == 0 ? D : zeros(D)
+
 ==(Da::Diagonal, Db::Diagonal) = Da.diag == Db.diag
 
 +(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag + Db.diag)

base/linalg/factorization.jl

@@ -191,17 +191,20 @@ convert{T}(::Type{QRPivoted{T}},A::QRPivoted) = QRPivoted(convert(AbstractMatrix
 convert{T}(::Type{Factorization{T}}, A::QRPivoted) = convert(QRPivoted{T}, A)
 
 function getindex(A::QR, d::Symbol)
-    d == :R && return triu(A.factors[1:minimum(size(A)),:])
+    m, n = size(A)
+    d == :R && return triu!(A.factors[1:min(m,n), 1:n])
     d == :Q && return QRPackedQ(A.factors,A.τ)
     throw(KeyError(d))
 end
 function getindex(A::QRCompactWY, d::Symbol)
-    d == :R && return triu(A.factors[1:minimum(size(A)),:])
+    m, n = size(A)
+    d == :R && return triu!(A.factors[1:min(m,n), 1:n])
     d == :Q && return QRCompactWYQ(A.factors,A.T)
     throw(KeyError(d))
 end
 function getindex{T}(A::QRPivoted{T}, d::Symbol)
-    d == :R && return triu(A.factors[1:minimum(size(A)),:])
+    m, n = size(A)
+    d == :R && return triu!(A.factors[1:min(m,n), 1:n])
     d == :Q && return QRPackedQ(A.factors,A.τ)
     d == :p && return A.jpvt
     if d == :P

base/linalg/generic.jl

@@ -22,15 +22,11 @@ scale!(s::Number, X::AbstractArray) = generic_scale!(X, s)
 
 cross(a::AbstractVector, b::AbstractVector) = [a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]]
 
-triu(M::AbstractMatrix) = triu(M,0)
-tril(M::AbstractMatrix) = tril(M,0)
-#triu{T}(M::AbstractMatrix{T}, k::Integer)
-#tril{T}(M::AbstractMatrix{T}, k::Integer)
+triu(M::AbstractMatrix) = triu!(copy(M))
+tril(M::AbstractMatrix) = tril!(copy(M))
 triu!(M::AbstractMatrix) = triu!(M,0)
 tril!(M::AbstractMatrix) = tril!(M,0)
 
-#diff(a::AbstractVector)
-#diff(a::AbstractMatrix, dim::Integer)
 diff(a::AbstractMatrix) = diff(a, 1)
 diff(a::AbstractVector) = [ a[i+1] - a[i] for i=1:length(a)-1 ]
 
@@ -45,10 +41,8 @@ end
 
 gradient(F::AbstractVector) = gradient(F, [1:length(F)])
 gradient(F::AbstractVector, h::Real) = gradient(F, [h*(1:length(F))])
-#gradient(F::AbstractVector, h::AbstractVector)
 
 diag(A::AbstractVector) = error("use diagm instead of diag to construct a diagonal matrix")
-#diag(A::AbstractMatrix)
 
 #diagm{T}(v::AbstractVecOrMat{T})
 

base/linalg/lu.jl

@@ -93,8 +93,12 @@ end
 
 function getindex{T,S<:StridedMatrix}(A::LU{T,S}, d::Symbol)
     m, n = size(A)
-    d == :L && return tril(A.factors[1:m, 1:min(m,n)], -1) + eye(T, m, min(m,n))
-    d == :U && return triu(A.factors[1:min(m,n),1:n])
+    if d == :L
+        L = tril!(A.factors[1:m, 1:min(m,n)])
+        for i = 1:min(m,n); L[i,i] = one(T); end
+        return L
+    end
+    d == :U && return triu!(A.factors[1:min(m,n), 1:n])
     d == :p && return ipiv2perm(A.ipiv, m)
     if d == :P
         p = A[:p]
@@ -144,7 +148,7 @@ end
 
 inv{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}) = @assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
 
-cond{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
+cond{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm((A[:L]*A[:U])[A[:p],:], p)))
 cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)
 
 # Tridiagonal

base/linalg/special.jl

@@ -4,7 +4,8 @@
 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::Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal))=Triangular(full(A))
+convert(::Type{Triangular}, A::Diagonal) = Triangular(full(A), :L)
+convert(::Type{Triangular}, A::Bidiagonal) = Triangular(full(A), A.isupper ? :U : :L)
 convert(::Type{Matrix}, D::Diagonal) = diagm(D.diag)
 
 function convert(::Type{Diagonal}, A::Union(Bidiagonal, SymTridiagonal))

test/linalg1.jl

@@ -94,7 +94,7 @@ debug && println("(Automatic) Square LU decomposition")
     lua   = factorize(a)
     l,u,p = lua[:L], lua[:U], lua[:p]
     @test_approx_eq l*u a[p,:]
-    @test_approx_eq l[invperm(p),:]*u a
+    @test_approx_eq (l*u)[invperm(p),:] a
     @test_approx_eq a * inv(lua) eye(n)
     @test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
 

test/linalg4.jl

@@ -31,8 +31,7 @@ for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
         b += im*convert(Matrix{elty}, randn(n, 2))
     end
 
-    for M in (triu(A), tril(A))
-        TM = Triangular(M)
+    for (M, TM) in ((triu(A), Triangular(A, :U)), (tril(A), Triangular(A, :L)))
 
         ##Idempotent tests #XXX - not implemented
         #for func in (conj, transpose, ctranspose)
@@ -65,8 +64,7 @@ for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
 
         #Binary operations
         B = convert(Matrix{elty}, randn(n, n))
-        for M2 in (triu(B), tril(B))
-            TM2 = Triangular(M2)
+        for (M2, TM2) in ((triu(B), Triangular(B, :U)), (tril(B), Triangular(B, :L)))
             for op in (*,) #+, - not implemented
                 @test_approx_eq full(op(TM, TM2)) op(M, M2)
                 @test_approx_eq full(op(TM, M2)) op(M, M2)
@@ -301,11 +299,11 @@ for newtype in [Tridiagonal, Matrix]
 end
 
 A=Tridiagonal(zeros(n-1), [1.0:n], zeros(n-1)) #morally Diagonal
-for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix]
+for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix]
     @test full(convert(newtype, A)) == full(A)
 end
 
-A=Triangular(full(Diagonal(a))) #morally Diagonal
+A=Triangular(full(Diagonal(a)), :L) #morally Diagonal
 for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix]
     @test full(convert(newtype, A)) == full(A)
 end