From 1c08f8df825f31fcb7a2df4093ae57c696cde26b Mon Sep 17 00:00:00 2001 From: kshyatt Date: Wed, 12 Aug 2015 20:58:02 -0700 Subject: [PATCH] Fixed symmetric --- base/linalg/symmetric.jl | 51 ++++++++++++++++++++++++++++++++++++---- test/linalg/symmetric.jl | 14 +++++++---- 2 files changed, 57 insertions(+), 8 deletions(-) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index bb37b51a5d06e..408d882e91613 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -57,10 +57,53 @@ ctranspose(A::Hermitian) = A trace(A::Hermitian) = real(trace(A.data)) #tril/triu -tril(A::Hermitian,k::Integer=0) = tril(A.data,k) -triu(A::Hermitian,k::Integer=0) = triu(A.data,k) -tril(A::Symmetric,k::Integer=0) = tril(A.data,k) -triu(A::Symmetric,k::Integer=0) = triu(A.data,k) +function tril(A::Hermitian, k::Integer=0) + if A.uplo == 'U' && k <= 0 + return tril!(A.data',k) + elseif A.uplo == 'U' && k > 0 + return tril!(A.data',-1) + tril!(triu(A.data),k) + elseif A.uplo == 'L' && k <= 0 + return tril(A.data,k) + else + return tril(A.data,-1) + tril!(triu!(A.data'),k) + end +end + +function tril(A::Symmetric, k::Integer=0) + if A.uplo == 'U' && k <= 0 + return tril!(A.data.',k) + elseif A.uplo == 'U' && k > 0 + return tril!(A.data.',-1) + tril!(triu(A.data),k) + elseif A.uplo == 'L' && k <= 0 + return tril(A.data,k) + else + return tril(A.data,-1) + tril!(triu!(A.data.'),k) + end +end + +function triu(A::Hermitian, k::Integer=0) + if A.uplo == 'U' && k >= 0 + return triu(A.data,k) + elseif A.uplo == 'U' && k < 0 + return triu(A.data,1) + triu!(tril!(A.data'),k) + elseif A.uplo == 'L' && k >= 0 + return triu!(A.data',k) + else + return triu!(A.data',1) + triu!(tril(A.data),k) + end +end + +function triu(A::Symmetric, k::Integer=0) + if A.uplo == 'U' && k >= 0 + return triu(A.data,k) + elseif A.uplo == 'U' && k < 0 + return triu(A.data,1) + triu!(tril!(A.data.'),k) + elseif A.uplo == 'L' && k >= 0 + return triu!(A.data.',k) + else + return triu!(A.data.',1) + triu!(tril(A.data),k) + end +end ## Matvec A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(y::StridedVector{T}, A::Symmetric{T,S}, x::StridedVector{T}) = BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y) diff --git a/test/linalg/symmetric.jl b/test/linalg/symmetric.jl index 7c8b3c6ec34d2..c1511ac42ef3f 100644 --- a/test/linalg/symmetric.jl +++ b/test/linalg/symmetric.jl @@ -60,10 +60,16 @@ let n=10 @test ctranspose(Hermitian(asym)) == asym #tril/triu - @test triu(Symmetric(asym),1) == triu(asym,1) - @test tril(Symmetric(asym),1) == tril(asym,1) - @test triu(Hermitian(asym),1) == triu(asym,1) - @test tril(Hermitian(asym),1) == tril(asym,1) + for di in -n:n + @test triu(Symmetric(a+a.'),di) == triu(a+a.',di) + @test tril(Symmetric(a+a.'),di) == tril(a+a.',di) + @test triu(Hermitian(asym),di) == triu(asym,di) + @test tril(Hermitian(asym),di) == tril(asym,di) + @test triu(Symmetric(a+a.',:L),di) == triu(a+a.',di) + @test tril(Symmetric(a+a.',:L),di) == tril(a+a.',di) + @test triu(Hermitian(asym,:L),di) == triu(asym,di) + @test tril(Hermitian(asym,:L),di) == tril(asym,di) + end eltya == BigFloat && continue # Revisit when implemented in julia d, v = eig(asym)