diff --git a/base/exports.jl b/base/exports.jl index 7e2f9340909f9..5a795de8d57c7 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -16,6 +16,7 @@ export Markdown, # Types + €, AbstractMatrix, AbstractSparseArray, AbstractSparseMatrix, diff --git a/base/linalg.jl b/base/linalg.jl index 4b34c00a804b9..1907f977593df 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -9,6 +9,7 @@ export BLAS, # Types + €, SymTridiagonal, Tridiagonal, Bidiagonal, @@ -193,6 +194,7 @@ const CHARU = 'U' const CHARL = 'L' char_uplo(uplo::Symbol) = uplo == :U ? CHARU : (uplo == :L ? CHARL : throw(ArgumentError("uplo argument must be either :U or :L"))) +immutable €{T} end include("linalg/exceptions.jl") include("linalg/generic.jl") diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index 8c71fca854567..f1afd0796be91 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -13,20 +13,20 @@ immutable CholeskyPivoted{T} <: Factorization{T} info::BlasInt end -function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U) - C, info = LAPACK.potrf!(char_uplo(uplo), A) - return @assertposdef Triangular{eltype(C),typeof(C),uplo,false}(C) info +function chol!{T<:BlasFloat,UpLo}(A::StridedMatrix{T}, ::Type{€{UpLo}}) + C, info = LAPACK.potrf!(char_uplo(UpLo), A) + return @assertposdef Triangular{eltype(C),typeof(C),UpLo,false}(C) info end -function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U) +function chol!{T,UpLo}(A::AbstractMatrix{T}, ::Type{€{UpLo}}) n = chksquare(A) @inbounds begin - if uplo == :L + if UpLo == :L for k = 1:n for i = 1:k - 1 A[k,k] -= A[k,i]*A[k,i]' end - A[k,k] = chol!(A[k,k], uplo) + A[k,k] = chol!(A[k,k], €{UpLo}) AkkInv = inv(A[k,k]') for j = 1:k for i = k + 1:n @@ -35,12 +35,12 @@ function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U) end end end - elseif uplo == :U + elseif UpLo == :U 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) + A[k,k] = chol!(A[k,k], €{UpLo}) AkkInv = inv(A[k,k]) for j = k + 1:n for i = 1:k - 1 @@ -50,21 +50,23 @@ function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U) end end else - throw(ArgumentError("uplo must be either :U or :L but was $(uplo)")) + throw(ArgumentError("UpLo must be either :U or :L but was $(UpLo)")) end end - return Triangular(A, uplo, false) + return Triangular(A, €{UpLo}, €{false}) end +chol!(A::AbstractMatrix) = chol!(A, €{:U}) -function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) - uplochar = char_uplo(uplo) +function cholfact!{T<:BlasFloat,UpLo}(A::StridedMatrix{T}, ::Type{€{UpLo}}; pivot=false, tol=0.0) + uplochar = char_uplo(UpLo) if pivot A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol) return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info) end - return Cholesky{T,typeof(A),uplo}(chol!(A, uplo).data) + return Cholesky{T,typeof(A),UpLo}(chol!(A, €{UpLo}).data) end -cholfact!(A::AbstractMatrix, uplo::Symbol=:U) = Cholesky{eltype(A),typeof(A),uplo}(chol!(A, uplo).data) +cholfact!{UpLo}(A::AbstractMatrix, ::Type{€{UpLo}}) = Cholesky{eltype(A),typeof(A),UpLo}(chol!(A, €{UpLo}).data) +cholfact!(A::AbstractMatrix) = cholfact!(A, €{:U}) function cholfact!{T<:BlasFloat,S,UpLo}(C::Cholesky{T,S,UpLo}) _, info = LAPACK.potrf!(char_uplo(UpLo), C.UL) @@ -72,27 +74,31 @@ function cholfact!{T<:BlasFloat,S,UpLo}(C::Cholesky{T,S,UpLo}) C end -cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol) -function cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) +cholfact{T<:BlasFloat,UpLo}(A::StridedMatrix{T}, ::Type{€{UpLo}}; pivot=false, tol=0.0) = cholfact!(copy(A), €{UpLo}, pivot=pivot, tol=tol) +function cholfact{T,UpLo}(A::StridedMatrix{T}, ::Type{€{UpLo}}; pivot=false, tol=0.0) S = promote_type(typeof(chol(one(T))),Float32) - S <: BlasFloat && return cholfact!(convert(AbstractMatrix{S}, A), uplo, pivot = pivot, tol = tol) + S <: BlasFloat && return cholfact!(convert(AbstractMatrix{S}, A), €{UpLo}, pivot = pivot, tol = tol) pivot && throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}")) - S != T && return cholfact!(convert(AbstractMatrix{S}, A), uplo) - return cholfact!(copy(A), uplo) + S != T && return cholfact!(convert(AbstractMatrix{S}, A), €{UpLo}) + return cholfact!(copy(A), €{UpLo}) end -function cholfact(x::Number, uplo::Symbol=:U) +function cholfact{UpLo}(x::Number, ::Type{€{UpLo}}) xf = fill(chol!(x), 1, 1) Cholesky{:U, eltype(xf), typeof(xf)}(xf) end +cholfact(A; pivot=false, tol=0.0) = cholfact(A, €{:U}; pivot = pivot, tol = tol) -chol(A::AbstractMatrix, uplo::Symbol=:U) = chol!(copy(A), uplo) -function chol!(x::Number, uplo::Symbol=:U) +chol{UpLo}(A::AbstractMatrix, ::Type{€{UpLo}}) = chol!(copy(A), €{UpLo}) +chol(A::AbstractMatrix) = chol!(copy(A), €{:U}) +function chol!(x::Number) rx = real(x) rx == abs(x) || throw(DomainError()) rxr = sqrt(rx) convert(promote_type(typeof(x), typeof(rxr)), rxr) end -chol(x::Number, uplo::Symbol=:U) = chol!(x, uplo) +chol!(x::Number, ::Type) = chol!(x) +chol(x::Number) = chol!(x) +chol(x::Number, ::Type) = chol!(x, €{:U}) function convert{Tnew,Told,S,UpLo}(::Type{Cholesky{Tnew}},C::Cholesky{Told,S,UpLo}) Cnew = convert(AbstractMatrix{Tnew}, C.UL) @@ -106,39 +112,33 @@ convert{T}(::Type{Factorization{T}}, C::Cholesky) = convert(Cholesky{T}, C) convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) = CholeskyPivoted(convert(AbstractMatrix{T},C.UL),C.uplo,C.piv,C.rank,C.tol,C.info) convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C) -full{T,S}(C::Cholesky{T,S,:U}) = C[:U]'C[:U] -full{T,S}(C::Cholesky{T,S,:L}) = C[:L]*C[:L]' +full{T,S}(C::Cholesky{T,S,:U}) = C[€{:U}]'C[€{:U}] +full{T,S}(C::Cholesky{T,S,:L}) = C[€{:L}]*C[€{:L}]' 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) - 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 == :p && return C.piv - if d == :P - n = size(C, 1) - P = zeros(T, n, n) - for i=1:n - P[C.piv[i],i] = one(T) - end - return P +getindex{T,S,UpLo}(C::Cholesky{T,S,UpLo}, ::Type{€{:U}}) = Triangular(UpLo == :U ? C.UL : C.UL', €{:U}) +getindex{T,S,UpLo}(C::Cholesky{T,S,UpLo}, ::Type{€{:L}}) = Triangular(UpLo == :L ? C.UL : C.UL', €{:L}) +getindex{T,S,UpLo}(C::Cholesky{T,S,UpLo}, ::Type{€{:UL}}) = Triangular(C.UL, €{UpLo}) +getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, ::Type{€{:U}}) = Triangular(C.uplo == 'U' ? C.UL : C.UL', €{:U}) +getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, ::Type{€{:L}}) = Triangular(C.uplo == 'L' ? C.UL : C.UL', €{:L}) +getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, ::Type{€{:p}}) = C.piv +function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, ::Type{€{:P}}) + n = size(C, 1) + P = zeros(T, n, n) + for i = 1:n + P[C.piv[i],i] = one(T) end - throw(KeyError(d)) + return P end show{T,S<:AbstractMatrix,UpLo}(io::IO, C::Cholesky{T,S,UpLo}) = (println("$(typeof(C)) with factor:");show(io,C[UpLo])) 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!(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)) function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T}) chkfullrank(C) @@ -156,8 +156,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!(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),:] 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 474a494c72005..96df25b472c44 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -307,16 +307,16 @@ 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))) - retmat = SchurF[:vectors]*R*SchurF[:vectors]' + 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], :U, false))) - SchurF[:vectors]*R*SchurF[:vectors]' + 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) sqrtm(a::Complex) = sqrt(a) @@ -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(Triangular(A, €{:U}, €{false})) elseif istril(Ac) - Ai = inv(Triangular(A, :L, false)) + Ai = inv(Triangular(A, €{:L}, €{false})) 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 Triangular(A, €{:L}) 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 Triangular(A, €{:U}) 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) : \(Triangular(A, €{:L}),B) end - istriu(A) && return \(Triangular(A, :U),B) + istriu(A) && return \(Triangular(A, €{:U}),B) return \(lufact(A),B) end return qrfact(A,pivot=eltype(A)<:BlasFloat)\B @@ -440,12 +440,12 @@ function pinv{T}(A::StridedMatrix{T}, tol::Real) end end SVD = svdfact(A, thin=true) - S = eltype(SVD[:S]) - Sinv = zeros(S, length(SVD[:S])) - index = SVD[:S] .> tol*maximum(SVD[:S]) - Sinv[index] = one(S) ./ SVD[:S][index] + S = eltype(SVD[€{:S}]) + Sinv = zeros(S, length(SVD[€{:S}])) + index = SVD[€{:S}] .> tol*maximum(SVD[€{:S}]) + Sinv[index] = one(S) ./ SVD[€{:S}][index] Sinv[find(!isfinite(Sinv))] = zero(S) - return SVD[:Vt]'scale(Sinv, SVD[:U]') + return SVD[€{:Vt}]'scale(Sinv, SVD[€{:U}]') end function pinv{T}(A::StridedMatrix{T}) tol = eps(real(float(one(T))))*maximum(size(A)) @@ -459,8 +459,8 @@ function null{T}(A::StridedMatrix{T}) m, n = size(A) (m == 0 || n == 0) && return eye(T, n) SVD = svdfact(A, thin=false) - indstart = sum(SVD[:S] .> max(m,n)*maximum(SVD[:S])*eps(eltype(SVD[:S]))) + 1 - return SVD[:V][:,indstart:end] + indstart = sum(SVD[€{:S}] .> max(m,n)*maximum(SVD[€{:S}])*eps(eltype(SVD[€{:S}]))) + 1 + return SVD[€{:V}][:,indstart:end] end null(a::StridedVector) = null(reshape(a, length(a), 1)) diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index c8531517afaa3..b559e6a21b9b0 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -64,9 +64,9 @@ qrfact(x::Number) = qrfact(fill(x,1,1)) function qr(A::Union(Number, AbstractMatrix); pivot::Bool=false, thin::Bool=true) F = qrfact(A, pivot=pivot) if pivot - full(F[:Q], thin=thin), F[:R], F[:p] + full(F[€{:Q}], thin=thin), F[€{:R}], F[€{:p}] else - full(F[:Q], thin=thin), F[:R] + full(F[€{:Q}], thin=thin), F[€{:R}] end end @@ -77,33 +77,31 @@ convert{T}(::Type{Factorization{T}}, A::QRCompactWY) = convert(QRCompactWY{T}, A convert{T}(::Type{QRPivoted{T}},A::QRPivoted) = QRPivoted(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ), A.jpvt) convert{T}(::Type{Factorization{T}}, A::QRPivoted) = convert(QRPivoted{T}, A) -function getindex(A::QR, d::Symbol) +function getindex(A::QR, ::Type{€{:R}}) 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)) + triu!(A.factors[1:min(m,n), 1:n]) end -function getindex(A::QRCompactWY, d::Symbol) +getindex(A::QR, ::Type{€{:Q}}) = QRPackedQ(A.factors, A.τ) +function getindex(A::QRCompactWY, ::Type{€{:R}}) 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)) + return triu!(A.factors[1:min(m,n), 1:n]) end -function getindex{T}(A::QRPivoted{T}, d::Symbol) +getindex(A::QRCompactWY, ::Type{€{:Q}}) = QRCompactWYQ(A.factors, A.T) +function getindex{T}(A::QRPivoted{T}, ::Type{€{:R}}) 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 - p = A[:p] - n = length(p) - P = zeros(T, n, n) - for i in 1:n - P[p[i],i] = one(T) - end - return P + return triu!(A.factors[1:min(m,n), 1:n]) +end +getindex{T}(A::QRPivoted{T}, ::Type{€{:Q}}) = QRPackedQ(A.factors, A.τ) +getindex{T}(A::QRPivoted{T}, ::Type{€{:p}}) = A.jpvt +function getindex{T}(A::QRPivoted{T}, ::Type{€{:P}}) + m, n = size(A) + p = A[€{:p}] + n = length(p) + P = zeros(T, n, n) + for i in 1:n + P[p[i],i] = one(T) end - throw(KeyError(d)) + return P end # Type-stable interface to get Q getq(A::QRCompactWY) = QRCompactWYQ(A.factors,A.T) @@ -265,8 +263,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))) -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))) +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) # Julia implementation similarly to xgelsy function A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) @@ -294,10 +292,11 @@ 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!(Triangular(C[1:rnk,1:rnk], €{:U}),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)) - return isa(A,QRPivoted) ? B[invperm(A[:p]::Vector{BlasInt}),:] : B[1:nA,:], rnk + B[1:nA,:] = sub(B, 1:nA, :)[invperm(A[€{:p}]::Vector{BlasInt}),:] + return B, rnk end A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedVector{T}) = vec(A_ldiv_B!(A,reshape(B,length(B),1))) A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedVecOrMat{T}) = A_ldiv_B!(A, B, maximum(size(A))*eps(real(float(one(eltype(B))))))[1] @@ -305,8 +304,8 @@ function A_ldiv_B!{T}(A::QR{T},B::StridedMatrix{T}) m, n = size(A) minmn = min(m,n) mB, nB = size(B) - Ac_mul_B!(A[:Q],sub(B,1:m,1:nB)) # Reconsider when arrayviews are merged. - R = A[:R] + Ac_mul_B!(A[€{:Q}],sub(B,1:m,1:nB)) # Reconsider when arrayviews are merged. + R = A[€{:R}] @inbounds begin if n > m # minimum norm solution τ = zeros(T,m) @@ -350,22 +349,32 @@ function A_ldiv_B!{T}(A::QR{T},B::StridedMatrix{T}) end end end - return B[1:n,:] + return B end A_ldiv_B!(A::QR, B::StridedVector) = A_ldiv_B!(A, reshape(B, length(B), 1))[:] -A_ldiv_B!(A::QRPivoted, B::StridedVector) = A_ldiv_B!(QR(A.factors,A.τ),B)[invperm(A.jpvt)] -A_ldiv_B!(A::QRPivoted, B::StridedMatrix) = A_ldiv_B!(QR(A.factors,A.τ),B)[invperm(A.jpvt),:] +function A_ldiv_B!(A::QRPivoted, b::StridedVector) + A_ldiv_B!(QR(A.factors,A.τ), b) + b[1:size(A.factors, 2)] = sub(b, 1:size(A.factors, 2))[invperm(A.jpvt)] + b +end +function A_ldiv_B!(A::QRPivoted, B::StridedMatrix) + A_ldiv_B!(QR(A.factors, A.τ), B) + B[1:size(A.factors, 2),:] = sub(B, 1:size(A.factors, 2), :)[invperm(A.jpvt)] + B +end function \{TA,Tb}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),b::StridedVector{Tb}) S = promote_type(TA,Tb) m,n = size(A) m == length(b) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has length $(length(b))")) - n > m ? A_ldiv_B!(convert(Factorization{S},A),[b,zeros(S,n-m)]) : A_ldiv_B!(convert(Factorization{S},A), S == Tb ? copy(b) : convert(AbstractVector{S}, b)) + x = n > m ? A_ldiv_B!(convert(Factorization{S},A),[b,zeros(S,n-m)]) : A_ldiv_B!(convert(Factorization{S},A), S == Tb ? copy(b) : convert(AbstractVector{S}, b)) + return length(x) > n ? x[1:n] : x end function \{TA,TB}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),B::StridedMatrix{TB}) S = promote_type(TA,TB) m,n = size(A) m == size(B,1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(B,1)) rows")) - n > m ? A_ldiv_B!(convert(Factorization{S},A),[B;zeros(S,n-m,size(B,2))]) : A_ldiv_B!(convert(Factorization{S},A), S == TB ? copy(B) : convert(AbstractMatrix{S}, B)) + X = n > m ? A_ldiv_B!(convert(Factorization{S},A),[B;zeros(S,n-m,size(B,2))]) : A_ldiv_B!(convert(Factorization{S},A), S == TB ? copy(B) : convert(AbstractMatrix{S}, B)) + return size(X, 1) > n ? X[1:n,:] : X end ##TODO: Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly @@ -390,11 +399,8 @@ end HessenbergQ(A::Hessenberg) = HessenbergQ(A.factors, A.τ) size(A::HessenbergQ, args...) = size(A.factors, args...) -function getindex(A::Hessenberg, d::Symbol) - d == :Q && return HessenbergQ(A) - d == :H && return triu(A.factors, -1) - throw(KeyError(d)) -end +getindex(A::Hessenberg, ::Type{€{:Q}}) = HessenbergQ(A) +getindex(A::Hessenberg, ::Type{€{:H}}) = triu(A.factors, -1) full(A::HessenbergQ) = LAPACK.orghr!(1, size(A.factors, 1), copy(A.factors), A.τ) @@ -418,11 +424,8 @@ immutable GeneralizedEigen{T,V} <: Factorization{T} vectors::Matrix{T} end -function getindex(A::Union(Eigen,GeneralizedEigen), d::Symbol) - d == :values && return A.values - d == :vectors && return A.vectors - throw(KeyError(d)) -end +getindex(A::Union(Eigen,GeneralizedEigen), ::Type{€{:values}}) = A.values +getindex(A::Union(Eigen,GeneralizedEigen), ::Type{€{:vectors}}) = A.vectors isposdef(A::Union(Eigen,GeneralizedEigen)) = all(A.values .> 0) @@ -462,10 +465,10 @@ eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1)) # end function eig(A::Union(Number, AbstractMatrix), args...; kwargs...) F = eigfact(A, args..., kwargs...) - F[:values], F[:vectors] + F[€{:values}], F[€{:vectors}] end #Calculates eigenvectors -eigvecs(A::Union(Number, AbstractMatrix), args...; kwargs...) = eigfact(A, args...; kwargs...)[:vectors] +eigvecs(A::Union(Number, AbstractMatrix), args...; kwargs...) = eigfact(A, args...; kwargs...)[€{:vectors}] function eigvals!{T<:BlasReal}(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) issym(A) && return eigvals!(Symmetric(A)) @@ -558,13 +561,10 @@ function svd(A::Union(Number, AbstractArray); thin::Bool=true) F.U, F.S, F.Vt' end -function getindex(F::SVD, d::Symbol) - d == :U && return F.U - d == :S && return F.S - d == :Vt && return F.Vt - d == :V && return F.Vt' - throw(KeyError(d)) -end +getindex(F::SVD, d::Type{€{:U}}) = F.U +getindex(F::SVD, d::Type{€{:S}}) = F.S +getindex(F::SVD, d::Type{€{:Vt}}) = F.Vt +getindex(F::SVD, d::Type{€{:V}}) = F.Vt' svdvals!{T<:BlasFloat}(A::StridedMatrix{T}) = any([size(A)...].==0) ? zeros(T, 0) : LAPACK.gesdd!('N', A)[2] svdvals{T<:BlasFloat}(A::StridedMatrix{T}) = svdvals!(copy(A)) @@ -601,39 +601,39 @@ svdfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(F function svd(A::AbstractMatrix, B::AbstractMatrix) F = svdfact(A, B) - F[:U], F[:V], F[:Q], F[:D1], F[:D2], F[:R0] -end - -function getindex{T}(obj::GeneralizedSVD{T}, d::Symbol) - d == :U && return obj.U - d == :V && return obj.V - d == :Q && return obj.Q - (d == :alpha || d == :a) && return obj.a - (d == :beta || d == :b) && return obj.b - (d == :vals || d == :S) && return obj.a[1:obj.k + obj.l] ./ obj.b[1:obj.k + obj.l] - if d == :D1 - m = size(obj.U, 1) - if m - obj.k - obj.l >= 0 - return [eye(T, obj.k) zeros(T, obj.k, obj.l); zeros(T, obj.l, obj.k) diagm(obj.a[obj.k + 1:obj.k + obj.l]); zeros(T, m - obj.k - obj.l, obj.k + obj.l)] - else - return [eye(T, m, obj.k) [zeros(T, obj.k, m - obj.k); diagm(obj.a[obj.k + 1:m])] zeros(T, m, obj.k + obj.l - m)] - end - end - if d == :D2 - m = size(obj.U, 1) - p = size(obj.V, 1) - if m - obj.k - obj.l >= 0 - return [zeros(T, obj.l, obj.k) diagm(obj.b[obj.k + 1:obj.k + obj.l]); zeros(T, p - obj.l, obj.k + obj.l)] - else - return [zeros(T, p, obj.k) [diagm(obj.b[obj.k + 1:m]); zeros(T, obj.k + p - m, m - obj.k)] [zeros(T, m - obj.k, obj.k + obj.l - m); eye(T, obj.k + p - m, obj.k + obj.l - m)]] - end + F[€{:U}], F[€{:V}], F[€{:Q}], F[€{:D1}], F[€{:D2}], F[€{:R0}] +end + +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:U}}) = obj.U +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:V}}) = obj.V +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:Q}}) = obj.Q +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:α}}) = obj.a +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:a}}) = obj.a +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:β}}) = obj.b +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:b}}) = obj.b +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:vals}}) = obj.a[1:obj.k + obj.l] ./ obj.b[1:obj.k + obj.l] +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:S}}) = obj.a[1:obj.k + obj.l] ./ obj.b[1:obj.k + obj.l] +function getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:D1}}) + m = size(obj.U, 1) + if m - obj.k - obj.l >= 0 + return [eye(T, obj.k) zeros(T, obj.k, obj.l); zeros(T, obj.l, obj.k) diagm(obj.a[obj.k + 1:obj.k + obj.l]); zeros(T, m - obj.k - obj.l, obj.k + obj.l)] + else + return [eye(T, m, obj.k) [zeros(T, obj.k, m - obj.k); diagm(obj.a[obj.k + 1:m])] zeros(T, m, obj.k + obj.l - m)] end - d == :R && return obj.R - if d == :R0 - n = size(obj.Q, 1) - return [zeros(T, obj.k + obj.l, n - obj.k - obj.l) obj.R] +end +function getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:D2}}) + m = size(obj.U, 1) + p = size(obj.V, 1) + if m - obj.k - obj.l >= 0 + return [zeros(T, obj.l, obj.k) diagm(obj.b[obj.k + 1:obj.k + obj.l]); zeros(T, p - obj.l, obj.k + obj.l)] + else + return [zeros(T, p, obj.k) [diagm(obj.b[obj.k + 1:m]); zeros(T, obj.k + p - m, m - obj.k)] [zeros(T, m - obj.k, obj.k + obj.l - m); eye(T, obj.k + p - m, obj.k + obj.l - m)]] end - throw(KeyError(d)) +end +getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:R}}) = obj.R +function getindex{T}(obj::GeneralizedSVD{T}, ::Type{€{:R0}}) + n = size(obj.Q, 1) + return [zeros(T, obj.k + obj.l, n - obj.k - obj.l) obj.R] end function svdvals!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) @@ -653,21 +653,20 @@ schurfact!{T<:BlasFloat}(A::StridedMatrix{T}) = Schur(LinAlg.LAPACK.gees!('V', A schurfact{T<:BlasFloat}(A::StridedMatrix{T}) = schurfact!(copy(A)) schurfact{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? schurfact!(convert(AbstractMatrix{S},A)) : schurfact!(copy(A))) -function getindex(F::Schur, d::Symbol) - (d == :T || d == :Schur) && return F.T - (d == :Z || d == :vectors) && return F.Z - d == :values && return F.values - throw(KeyError(d)) -end +getindex(F::Schur, ::Type{€{:T}}) = F.T +getindex(F::Schur, ::Type{€{:Schur}}) = F.T +getindex(F::Schur, ::Type{€{:Z}}) = F.Z +getindex(F::Schur, ::Type{€{:vectors}}) = F.Z +getindex(F::Schur, ::Type{€{:values}}) = F.values function schur(A::AbstractMatrix) SchurF = schurfact(A) - SchurF[:T], SchurF[:Z], SchurF[:values] + SchurF[€{:T}], SchurF[€{:Z}], SchurF[€{:values}] end ordschur!{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = Schur(LinAlg.LAPACK.trsen!(select, T , Q)...) ordschur{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = ordschur!(copy(Q), copy(T), select) -ordschur!{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = (res=ordschur!(schur.Z, schur.T, select); schur[:values][:]=res[:values]; res) +ordschur!{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = (res=ordschur!(schur.Z, schur.T, select); schur[€{:values}][:]=res[€{:values}]; res) ordschur{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = ordschur(schur.Z, schur.T, select) immutable GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty} @@ -683,20 +682,20 @@ schurfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = Generalized schurfact{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = schurfact!(copy(A),copy(B)) schurfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); schurfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B))) -function getindex(F::GeneralizedSchur, d::Symbol) - d == :S && return F.S - d == :T && return F.T - d == :alpha && return F.alpha - d == :beta && return F.beta - d == :values && return F.alpha./F.beta - (d == :Q || d == :left) && return F.Q - (d == :Z || d == :right) && return F.Z - throw(KeyError(d)) -end +getindex(F::GeneralizedSchur, ::Type{€{:S}}) = F.S +getindex(F::GeneralizedSchur, ::Type{€{:T}}) = F.T +getindex(F::GeneralizedSchur, ::Type{€{:α}}) = F.alpha +getindex(F::GeneralizedSchur, ::Type{€{:β}}) = F.beta +getindex(F::GeneralizedSchur, ::Type{€{:values}}) = F.alpha./F.beta +getindex(F::GeneralizedSchur, ::Type{€{:Q}}) = F.Q +getindex(F::GeneralizedSchur, ::Type{€{:L}}) = F.Q +getindex(F::GeneralizedSchur, ::Type{€{:Z}}) = F.Z +getindex(F::GeneralizedSchur, ::Type{€{:R}}) = F.Z + function schur(A::AbstractMatrix, B::AbstractMatrix) SchurF = schurfact(A, B) - SchurF[:S], SchurF[:T], SchurF[:Q], SchurF[:Z] + SchurF[€{:S}], SchurF[€{:T}], SchurF[€{:Q}], SchurF[€{:Z}] end ### General promotion rules diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index 0764d439d8900..d305eda0afa29 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -424,7 +424,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(Triangular(A, €{:U}, €{false})) return det(lufact(A)) end det(x::Number) = x diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl index c07806b885efd..0dd1e1a0d77e2 100644 --- a/base/linalg/lu.jl +++ b/base/linalg/lu.jl @@ -69,7 +69,7 @@ lufact(F::LU) = F lu(x::Number) = (one(x), x, 1) function lu(A::AbstractMatrix; pivot = true) F = lufact(A, pivot = pivot) - F[:L], F[:U], F[:p] + F[€{:L}], F[€{:U}], F[€{:p}] end function convert{T}(::Type{LU{T}}, F::LU) @@ -91,29 +91,33 @@ function ipiv2perm{T}(v::AbstractVector{T}, maxi::Integer) return p end -function getindex{T,S<:StridedMatrix}(A::LU{T,S}, d::Symbol) +function getindex{T,S<:StridedMatrix}(A::LU{T,S}, ::Type{€{:L}}) m, n = size(A) - 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] - P = zeros(T, m, m) - for i in 1:m - P[i,p[i]] = one(T) - end - return P + 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 +function getindex{T,S<:StridedMatrix}(A::LU{T,S}, ::Type{€{:U}}) + m, n = size(A) + triu!(A.factors[1:min(m,n), 1:n]) +end +function getindex{T,S<:StridedMatrix}(A::LU{T,S}, ::Type{€{:p}}) + m, n = size(A) + ipiv2perm(A.ipiv, m) +end +function getindex{T,S<:StridedMatrix}(A::LU{T,S}, ::Type{€{:P}}) + m, n = size(A) + p = A[€{:p}] + P = zeros(T, m, m) + for i in 1:m + P[i,p[i]] = one(T) end - throw(KeyError(d)) + return P 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!(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)),:])) 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 @@ -148,8 +152,8 @@ 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[:U])[A[:p],:], p))) -cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),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 diff --git a/base/linalg/special.jl b/base/linalg/special.jl index 4113e49446d95..ebad8966643b4 100644 --- a/base/linalg/special.jl +++ b/base/linalg/special.jl @@ -4,8 +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::Diagonal) = Triangular(full(A), :L) -convert(::Type{Triangular}, A::Bidiagonal) = Triangular(full(A), A.isupper ? :U : :L) +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)) diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index e82d2ff666582..25d5acabedc81 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -2,15 +2,19 @@ immutable Triangular{T,S<:AbstractMatrix,UpLo,IsUnit} <: AbstractMatrix{T} data::S end -function Triangular{T}(A::AbstractMatrix{T}, uplo::Symbol, isunit::Bool=false) +function Triangular{T,UpLo,UnitDiagonal}(A::AbstractMatrix{T}, ::Type{€{UpLo}}, ::Type{€{UnitDiagonal}}) chksquare(A) - Triangular{T,typeof(A),uplo,isunit}(A) + Triangular{T,typeof(A),UpLo,UnitDiagonal}(A) +end +function Triangular{T,UpLo}(A::AbstractMatrix{T}, ::Type{€{UpLo}}) + chksquare(A) + Triangular{T,typeof(A),UpLo,false}(A) end size(A::Triangular, args...) = size(A.data, args...) 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{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}) B = Array(Tret, size(A, 1), size(A, 1)) @@ -81,21 +85,21 @@ function -{T, S, UpLo}(A::Triangular{T, S, UpLo, true}) 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, 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, 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) ###################### @@ -166,39 +170,39 @@ eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:L,true}) = (for i = 1: # Generic routines # #################### -(*){T,S,UpLo}(A::Triangular{T,S,UpLo,false}, x::Number) = Triangular(A.data*x, UpLo) +(*){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) B = A.data*x for i = 1:size(A, 1) B[i,i] = x end - Triangular(B, UpLo) + Triangular(B, €{UpLo}) end -(*){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,false}) = Triangular(x*A.data, UpLo) +(*){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}) B = x*A.data for i = 1:size(A, 1) B[i,i] = x end - Triangular(B, UpLo) + Triangular(B, €{UpLo}) end -(/){T,S,UpLo}(A::Triangular{T,S,UpLo,false}, x::Number) = Triangular(A.data/x, UpLo) +(/){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) B = A.data*x invx = inv(x) for i = 1:size(A, 1) B[i,i] = x end - Triangular(B, UpLo) + Triangular(B, €{UpLo}) end -(\){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,false}) = Triangular(x\A.data, UpLo) +(\){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}) B = x\A.data invx = inv(x) for i = 1:size(A, 1) B[i,i] = invx end - Triangular(B, UpLo) + Triangular(B, €{UpLo}) end ## Generic triangular multiplication @@ -504,7 +508,7 @@ function sqrtm{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) r==0 || (R[i,j] = r / (R[i,i] + R[j,j])) end end - return Triangular(R, :U, IsUnit) + return Triangular(R, €{:U}, €{IsUnit}) else # UpLo == :L #Not the usual case return sqrtm(A.').' end diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl index bbb8be43c1734..439b73c496b21 100644 --- a/test/linalg/triangular.jl +++ b/test/linalg/triangular.jl @@ -8,8 +8,8 @@ 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 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] # 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) @@ -18,8 +18,8 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) 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) + 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)) @@ -29,9 +29,9 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) # 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) + @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) + @test_throws InexactError convert(Triangular{elty2}, A1) == Triangular(convert(Matrix{elty2}, A1.data), €{uplo1}, €{isunit1}) end @test convert(Matrix, A1) == full(A1) @@ -39,7 +39,7 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) @test full!(copy(A1)) == full(A1) # fill! - @test full!(fill!(copy(A1), 1)) == full(Triangular(ones(size(A1)...), uplo1, isunit1)) + @test full!(fill!(copy(A1), 1)) == full(Triangular(ones(size(A1)...), €{uplo1}, €{isunit1})) # similar @test isa(similar(A1), Triangular{elty1, Matrix{elty1}, uplo1, isunit1}) @@ -170,14 +170,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 = Triangular(lufact(A)[€{:U}], €{: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(Triangular(A, €{:U}), x, b)[1][i] end end debug && println("Test forward error [JIN 5705] if this is not a BigFloat") @@ -198,14 +198,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 = Triangular(lufact(A)[€{:U}], €{: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(Triangular(A, €{:U}), x, b)[1][i] end end diff --git a/test/linalg1.jl b/test/linalg1.jl index 31057b0d3ca93..8aafab2446f43 100644 --- a/test/linalg1.jl +++ b/test/linalg1.jl @@ -39,7 +39,7 @@ debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") debug && println("(Automatic) upper Cholesky factor") capd = factorize(apd) - r = capd[:U] + r = capd[€{:U}] κ = cond(apd, 1) #condition number #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 @@ -64,9 +64,9 @@ debug && println("(Automatic) upper Cholesky factor") @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow debug && println("lower Cholesky factor") - lapd = cholfact(apd, :L) + lapd = cholfact(apd, €{:L}) @test_approx_eq full(lapd) apd - l = lapd[:L] + l = lapd[€{:L}] @test_approx_eq l*l' apd debug && println("pivoted Choleksy decomposition") @@ -95,7 +95,7 @@ debug && println("Bunch-Kaufman factors of a pos-def matrix") debug && println("(Automatic) Square LU decomposition") κ = cond(a,1) lua = factorize(a) - l,u,p = lua[:L], lua[:U], lua[:p] + l,u,p = lua[€{:L}], lua[€{:U}], lua[€{:p}] @test_approx_eq l*u a[p,:] @test_approx_eq (l*u)[invperm(p),:] a @test_approx_eq a * inv(lua) eye(n) @@ -103,15 +103,15 @@ debug && println("(Automatic) Square LU decomposition") debug && println("Thin LU") lua = lufact(a[:,1:5]) - @test_approx_eq lua[:L]*lua[:U] lua[:P]*a[:,1:5] + @test_approx_eq lua[€{:L}]*lua[€{:U}] lua[€{:P}]*a[:,1:5] debug && println("Fat LU") lua = lufact(a[1:5,:]) - @test_approx_eq lua[:L]*lua[:U] lua[:P]*a[1:5,:] + @test_approx_eq lua[€{:L}]*lua[€{:U}] lua[€{:P}]*a[1:5,:] debug && println("QR decomposition (without pivoting)") qra = qrfact(a, pivot=false) - q,r = qra[:Q], qra[:R] + q,r = qra[€{:Q}], qra[€{:R}] @test_approx_eq q'*full(q, thin=false) eye(n) @test_approx_eq q*full(q, thin=false)' eye(n) @test_approx_eq q*r a @@ -119,8 +119,8 @@ debug && println("QR decomposition (without pivoting)") debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats qrpa = factorize(a[1:5,:]) - q,r = qrpa[:Q], qrpa[:R] - if isa(qrpa,QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia + q,r = qrpa[€{:Q}], qrpa[€{:R}] + if isa(qrpa,QRPivoted) p = qrpa[€{:p}] end # Reconsider if pivoted QR gets implemented in julia @test_approx_eq q'*full(q, thin=false) eye(5) @test_approx_eq q*full(q, thin=false)' eye(5) @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:5,p] : a[1:5,:] @@ -129,8 +129,8 @@ debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is onl debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats qrpa = factorize(a[:,1:5]) - q,r = qrpa[:Q], qrpa[:R] - if isa(qrpa, QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia + q,r = qrpa[€{:Q}], qrpa[€{:R}] + if isa(qrpa, QRPivoted) p = qrpa[€{:p}] end # Reconsider if pivoted QR gets implemented in julia @test_approx_eq q'*full(q, thin=false) eye(n) @test_approx_eq q*full(q, thin=false)' eye(n) @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:5] @@ -142,8 +142,8 @@ debug && println("symmetric eigen-decomposition") @test_approx_eq asym*v[:,1] d[1]*v[:,1] @test_approx_eq v*Diagonal(d)*v' asym @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1])) - @test_approx_eq abs(eigfact(Hermitian(asym), 1:2)[:vectors]'v[:,1:2]) eye(eltya, 2) - @test_approx_eq abs(eigfact(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2]))[:vectors]'v[:,1:2]) eye(eltya, 2) + @test_approx_eq abs(eigfact(Hermitian(asym), 1:2)[€{:vectors}]'v[:,1:2]) eye(eltya, 2) + @test_approx_eq abs(eigfact(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2]))[€{:vectors}]'v[:,1:2]) eye(eltya, 2) @test_approx_eq eigvals(Hermitian(asym), 1:2) d[1:2] @test_approx_eq eigvals(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) d[1:2] end @@ -158,26 +158,26 @@ debug && println("symmetric generalized eigenproblem") if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia a610 = a[:,6:10] f = eigfact(asym[1:5,1:5], a610'a610) - @test_approx_eq asym[1:5,1:5]*f[:vectors] scale(a610'a610*f[:vectors], f[:values]) - @test_approx_eq f[:values] eigvals(asym[1:5,1:5], a610'a610) - @test_approx_eq_eps prod(f[:values]) prod(eigvals(asym[1:5,1:5]/(a610'a610))) 200ε + @test_approx_eq asym[1:5,1:5]*f[€{:vectors}] scale(a610'a610*f[€{:vectors}], f[€{:values}]) + @test_approx_eq f[€{:values}] eigvals(asym[1:5,1:5], a610'a610) + @test_approx_eq_eps prod(f[€{:values}]) prod(eigvals(asym[1:5,1:5]/(a610'a610))) 200ε end debug && println("Non-symmetric generalized eigenproblem") if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia f = eigfact(a[1:5,1:5], a[6:10,6:10]) - @test_approx_eq a[1:5,1:5]*f[:vectors] scale(a[6:10,6:10]*f[:vectors], f[:values]) - @test_approx_eq f[:values] eigvals(a[1:5,1:5], a[6:10,6:10]) - @test_approx_eq_eps prod(f[:values]) prod(eigvals(a[1:5,1:5]/a[6:10,6:10])) 50000ε + @test_approx_eq a[1:5,1:5]*f[€{:vectors}] scale(a[6:10,6:10]*f[€{:vectors}], f[€{:values}]) + @test_approx_eq f[€{:values}] eigvals(a[1:5,1:5], a[6:10,6:10]) + @test_approx_eq_eps prod(f[€{:values}]) prod(eigvals(a[1:5,1:5]/a[6:10,6:10])) 50000ε end debug && println("Schur") if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia f = schurfact(a) - @test_approx_eq f[:vectors]*f[:Schur]*f[:vectors]' a - @test_approx_eq sort(real(f[:values])) sort(real(d)) - @test_approx_eq sort(imag(f[:values])) sort(imag(d)) - @test istriu(f[:Schur]) || iseltype(a,Real) + @test_approx_eq f[€{:vectors}]*f[€{:Schur}]*f[€{:vectors}]' a + @test_approx_eq sort(real(f[€{:values}])) sort(real(d)) + @test_approx_eq sort(imag(f[€{:values}])) sort(imag(d)) + @test istriu(f[€{:Schur}]) || iseltype(a,Real) end debug && println("Reorder Schur") @@ -188,30 +188,30 @@ debug && println("Reorder Schur") S = schurfact(ordschura) select = rand(range(0,2), n) O = ordschur(S, select) - bool(sum(select)) && @test_approx_eq S[:values][find(select)] O[:values][1:sum(select)] - @test_approx_eq O[:vectors]*O[:Schur]*O[:vectors]' ordschura + bool(sum(select)) && @test_approx_eq S[€{:values}][find(select)] O[€{:values}][1:sum(select)] + @test_approx_eq O[€{:vectors}]*O[€{:Schur}]*O[€{:vectors}]' ordschura end debug && println("Generalized Schur") if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia f = schurfact(a[1:5,1:5], a[6:10,6:10]) - @test_approx_eq f[:Q]*f[:S]*f[:Z]' a[1:5,1:5] - @test_approx_eq f[:Q]*f[:T]*f[:Z]' a[6:10,6:10] - @test istriu(f[:S]) || iseltype(a,Real) - @test istriu(f[:T]) || iseltype(a,Real) + @test_approx_eq f[€{:Q}]*f[€{:S}]*f[€{:Z}]' a[1:5,1:5] + @test_approx_eq f[€{:Q}]*f[€{:T}]*f[€{:Z}]' a[6:10,6:10] + @test istriu(f[€{:S}]) || iseltype(a,Real) + @test istriu(f[€{:T}]) || iseltype(a,Real) end debug && println("singular value decomposition") if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia usv = svdfact(a) - @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a + @test_approx_eq usv[€{:U}]*scale(usv[€{:S}],usv[€{:Vt}]) a end debug && println("Generalized svd") if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia gsvd = svdfact(a,a[1:5,:]) - @test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a - @test_approx_eq gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' a[1:5,:] + @test_approx_eq gsvd[€{:U}]*gsvd[€{:D1}]*gsvd[€{:R}]*gsvd[€{:Q}]' a + @test_approx_eq gsvd[€{:V}]*gsvd[€{:D2}]*gsvd[€{:R}]*gsvd[€{:Q}]' a[1:5,:] end debug && println("Solve square general system of equations") diff --git a/test/linalg2.jl b/test/linalg2.jl index f9cff896fa847..0244d20501af3 100644 --- a/test/linalg2.jl +++ b/test/linalg2.jl @@ -187,7 +187,7 @@ for elty in (Float32, Float64, Complex64, Complex128) A_mul_B!(G, R) end end - @test_approx_eq abs(A) abs(hessfact(Ac)[:H]) + @test_approx_eq abs(A) abs(hessfact(Ac)[€{:H}]) @test_approx_eq norm(R*eye(elty, 10)) one(elty) end @@ -232,7 +232,7 @@ end a = convert(Matrix{Rational{BigInt}}, rand(1:10//1,n,n))/n b = rand(1:10,n,2) lua = factorize(a) -l,u,p = lua[:L], lua[:U], lua[:p] +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 a * inv(lua) eye(n) @@ -255,7 +255,7 @@ for elty in (Float32, Float64, Complex64, Complex128) -eps(real(one(elty)))/4 eps(real(one(elty)))/2 -1.0 0; -0.5 -0.5 0.1 1.0]) F = eigfact(A,permute=false,scale=false) - @test_approx_eq F[:vectors]*Diagonal(F[:values])/F[:vectors] A + @test_approx_eq F[€{:vectors}]*Diagonal(F[€{:values}])/F[€{:vectors}] A F = eigfact(A) # @test norm(F[:vectors]*Diagonal(F[:values])/F[:vectors] - A) > 0.01 end @@ -422,7 +422,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 = Triangular(randn(3,3), €{:L}) @test T\I == inv(T) @test I\A == A @test A\I == inv(A) @@ -484,7 +484,7 @@ end #Issue 7304 let A=[-√.5 -√.5; -√.5 √.5] - Q=full(qrfact(A)[:Q]) + Q=full(qrfact(A)[€{:Q}]) @test vecnorm(A-Q) < eps() end diff --git a/test/linalg3.jl b/test/linalg3.jl index d834377f2ff8a..9c17c1a9141b4 100644 --- a/test/linalg3.jl +++ b/test/linalg3.jl @@ -169,7 +169,7 @@ for elty in (Float32, Float64, Complex64, Complex128) @test_approx_eq expm(A5) eA5 # Hessenberg - @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty}, + @test_approx_eq hessfact(A1)[€{:H}] convert(Matrix{elty}, [4.000000000000000 -1.414213562373094 -1.414213562373095 -1.414213562373095 4.999999999999996 -0.000000000000000 0 -0.000000000000002 3.000000000000000]) diff --git a/test/linalg4.jl b/test/linalg4.jl index c13e3693d55a3..69700046c489c 100644 --- a/test/linalg4.jl +++ b/test/linalg4.jl @@ -173,17 +173,17 @@ for isupper in (true, false) 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 +A = Triangular(full(Diagonal(a)), €{:L}) #morally Diagonal for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] @test full(convert(newtype, A)) == full(A) end @@ -196,8 +196,8 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) A = randn(5,5) end A = convert(Matrix{elty}, A'A) - for ul in (:U, :L) - @test_approx_eq full(cholfact(A, ul)[ul]) full(invoke(Base.LinAlg.chol!, (AbstractMatrix,Symbol),copy(A), ul)) + for ul in (€{:U}, €{:L}) + @test_approx_eq full(cholfact(A, ul)[ul]) full(invoke(Base.LinAlg.chol!, (AbstractMatrix,Type{ul}),copy(A), ul)) end end