Fixed symmetric

JuliaLang · Aug 14, 2015 · 1c08f8d · 1c08f8d
1 parent 095e0c4
commit 1c08f8d
Show file tree

Hide file tree

Showing 2 changed files with 57 additions and 8 deletions.
diff --git 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
@@ -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)