Add scale and copy for triangular matrices

JuliaLang · Aug 20, 2015 · 7210458 · 7210458
1 parent b6654a7
commit 7210458
Show file tree

Hide file tree

Showing 2 changed files with 91 additions and 1 deletion.
diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl
@@ -33,7 +33,7 @@ for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTr
             return $t(B)
         end
 
-        copy{T,S}(A::$t{T,S}) = $t{T,S}(copy(A.data))
+        # copy{T,S}(A::$t{T,S}) = $t{T,S}(copy(A.data))
 
         big(A::$t) = $t(big(A.data))
 
@@ -283,6 +283,90 @@ function -(A::UnitUpperTriangular)
     UpperTriangular(Anew)
 end
 
+# copy and scale
+function copy!{T<:Union{UpperTriangular, UnitUpperTriangular}}(A::T, B::T)
+    n = size(B,1)
+    for j = 1:n
+        for i = 1:(isa(B, UnitUpperTriangular)?j-1:j)
+            @inbounds A[i,j] = B[i,j]
+        end
+    end
+end
+function copy!{T<:Union{LowerTriangular, UnitLowerTriangular}}(A::T, B::T)
+    n = size(B,1)
+    for j = 1:n
+        for i = (isa(B, UnitLowerTriangular)?j+1:j):n
+            @inbounds A[i,j] = B[i,j]
+        end
+    end
+end
+function copy{T<:AbstractTriangular}(A::T)
+    D = similar(A)
+    copy!(D, A)
+    return D
+end
+for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular)
+    @eval begin
+        function copy{T<:$(t){BigFloat}}(A::T)
+            D = $(t)(zeros(BigFloat, size(A)...))
+            copy!(D, A)
+            return D
+        end
+    end
+end
+
+function scale!{T<:Union{UpperTriangular, UnitUpperTriangular}}(A::UpperTriangular, B::T, c::Number)
+    n = chksquare(B)
+    for j = 1:n
+        if isa(B, UnitUpperTriangular)
+            @inbounds A[j,j] = c
+        end
+        for i = 1:(isa(B, UnitUpperTriangular)?j-1:j)
+            @inbounds A[i,j] = c * B[i,j]
+        end
+    end
+end
+function scale!{T<:Union{LowerTriangular, UnitLowerTriangular}}(A::LowerTriangular, B::T, c::Number)
+    n = chksquare(B)
+    for j = 1:n
+        if isa(B, UnitLowerTriangular)
+            @inbounds A[j,j] = c
+        end
+        for i = (isa(B, UnitLowerTriangular)?j+1:j):n
+            @inbounds A[i,j] = c * B[i,j]
+        end
+    end
+end
+for (t1, t2) in ([:LowerTriangular, :LowerTriangular],
+                 [:UnitLowerTriangular, :LowerTriangular],
+                 [:UpperTriangular,:UpperTriangular],
+                 [:UnitUpperTriangular,UpperTriangular])
+    @eval begin
+        function scale{R<:Real,C<:Complex}(A::$(t1){R}, c::C)
+            D = $(t2)(Array(promote_type(R,C), size(A)...))
+            scale!(D,A,c)
+            return D
+        end
+        function scale{T,S<:Number}(A::$(t1){T}, c::S)
+            D = $(t2)(Array(promote_type(T,S), size(A)...))
+            scale!(D,A,c)
+            return D
+        end
+
+
+        function scale{R<:BigFloat,C<:Complex}(A::$(t1){R}, c::C)
+            D = $(t2)(zeros(promote_type(R,C),size(A)...))
+            scale!(D,A,c)
+            return D
+        end
+        function scale{T<:BigFloat,S<:Number}(A::$(t1){T}, c::S)
+            D = $(t2)(zeros(promote_type(T,S),size(A)...))
+            scale!(D,A,c)
+            return D
+        end
+    end
+end
+
 # Binary operations
 +(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data + B.data)
 +(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(A.data + B.data)

diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl
@@ -125,6 +125,12 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
         @test full(-A1) == -full(A1)
 
         # Binary operations
+        @test copy(A1) == A1
+        @test copy(A1.') == A1.'
+        @test scale(A1,0.5) == 0.5*A1
+        @test scale(A1,0.5im) == 0.5im*A1
+        @test scale(A1.',0.5) == 0.5*A1.'
+        @test scale(A1.',0.5im) == 0.5im*A1.'
         @test A1*0.5 == full(A1)*0.5
         @test 0.5*A1 == 0.5*full(A1)
         @test A1/0.5 == full(A1)/0.5