Add + and -(Triangular, UniformScaling) methods and move UniformScaling

tests to separate file

base/linalg/uniformscaling.jl

@@ -24,14 +24,63 @@ one{T}(J::UniformScaling{T}) = one(UniformScaling{T})
 zero{T}(::Type{UniformScaling{T}}) = UniformScaling(zero(T))
 zero{T}(J::UniformScaling{T}) = zero(UniformScaling{T})
 
-+(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ+J2.λ)
-+{T}(B::BitArray{2},J::UniformScaling{T}) = bitunpack(B) + J
-+(J::UniformScaling, B::BitArray{2}) = J + bitunpack(B)
-+(J::UniformScaling, A::AbstractMatrix) = A + J
-+(J::UniformScaling, x::Number) = J.λ + x
-+(x::Number, J::UniformScaling) = x + J.λ
-
-function +{TA,TJ}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
+(+)(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ+J2.λ)
+(+){T}(B::BitArray{2},J::UniformScaling{T}) = bitunpack(B) + J
+(+)(J::UniformScaling, B::BitArray{2}) = J + bitunpack(B)
+(+)(J::UniformScaling, A::AbstractMatrix) = A + J
+(+)(J::UniformScaling, x::Number) = J.λ + x
+(+)(x::Number, J::UniformScaling) = x + J.λ
+
+(-)(J::UniformScaling) = UniformScaling(-J.λ)
+(-)(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ-J2.λ)
+(-)(B::BitArray{2}, J::UniformScaling) = bitunpack(B) - J
+(-)(J::UniformScaling, B::BitArray{2}) = J - bitunpack(B)
+(-)(J::UniformScaling, x::Number) = J.λ - x
+(-)(x::Number, J::UniformScaling) = x - J.λ
+
+for (t1, t2) in ((:UnitUpperTriangular, :UpperTriangular),
+ (:UnitLowerTriangular, :LowerTriangular))
+ for op in (:+,:-)
+ @eval begin
+ ($op)(UL::$t2, J::UniformScaling) = ($t2)(($op)(UL.data, J))
+
+ function ($op)(UL::$t1, J::UniformScaling)
+ ULnew = copy_oftype(UL.data, promote_type(eltype(UL), eltype(J)))
+ for i = 1:size(ULnew, 1)
+ ULnew[i,i] = ($op)(1, J.λ)
+ end
+ return ($t2)(ULnew)
+ end
+ end
+ end
+
+ function (-)(J::UniformScaling, UL::Union{UpperTriangular,UnitUpperTriangular})
+ ULnew = similar(UL, promote_type(eltype(J), eltype(UL)))
+ n = size(ULnew, 1)
+ ULold = UL.data
+ for j = 1:n
+ for i = 1:j - 1
+ ULnew[i,j] = -ULold[i,j]
+ end
+ ULnew[j,j] = J.λ - ifelse(isa(UL, UnitUpperTriangular), 1, ULold[j,j])
+ end
+ return UpperTriangular(ULnew)
+ end
+ function (-)(J::UniformScaling, UL::Union{LowerTriangular,UnitLowerTriangular})
+ ULnew = similar(UL, promote_type(eltype(J), eltype(UL)))
+ n = size(ULnew, 1)
+ ULold = UL.data
+ for j = 1:n
+ ULnew[j,j] = J.λ - ifelse(isa(UL, UnitLowerTriangular), 1, ULold[j,j])
+ for i = j + 1:n
+ ULnew[i,j] = -ULold[i,j]
+ end
+ end
+ return LowerTriangular(ULNew)
+ end
+end
+
+function (+){TA,TJ}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
  n = chksquare(A)
  B = similar(A, promote_type(TA,TJ))
  copy!(B,A)
@@ -41,14 +90,7 @@ function +{TA,TJ}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
  B
 end
 
--(J::UniformScaling) = UniformScaling(-J.λ)
--(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ-J2.λ)
--(B::BitArray{2}, J::UniformScaling) = bitunpack(B) - J
--(J::UniformScaling, B::BitArray{2}) = J - bitunpack(B)
--(J::UniformScaling, x::Number) = J.λ - x
--(x::Number, J::UniformScaling) = x - J.λ
-
-function -{TA,TJ<:Number}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
+function (-){TA,TJ<:Number}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
  n = chksquare(A)
  B = similar(A, promote_type(TA,TJ))
  copy!(B, A)
@@ -57,11 +99,11 @@ function -{TA,TJ<:Number}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
  end
  B
 end
-function -{TA,TJ<:Number}(J::UniformScaling{TJ}, A::AbstractMatrix{TA})
+function (-){TA,TJ<:Number}(J::UniformScaling{TJ}, A::AbstractMatrix{TA})
  n = chksquare(A)
- B = -A
- @inbounds for i = 1:n
- B[i,i] += J.λ
+ B = convert(AbstractMatrix{promote_type(TJ,TA)}, -A)
+ @inbounds for j = 1:n
+ B[j,j] += J.λ
  end
  B
 end

test/choosetests.jl

@@ -65,7 +65,8 @@ function choosetests(choices = [])
  "linalg/lapack", "linalg/triangular", "linalg/tridiag",
  "linalg/bidiag", "linalg/diagonal", "linalg/dense",
  "linalg/pinv", "linalg/givens", "linalg/cholesky",
- "linalg/lu", "linalg/symmetric", "linalg/generic"]
+ "linalg/lu", "linalg/symmetric", "linalg/generic",
+ "linalg/uniformscaling"]
  if Base.USE_GPL_LIBS
  push!(linalgtests, "linalg/arnoldi")
  end

test/libdl.jl

@@ -1,3 +1,5 @@
+# This file is a part of Julia. License is MIT: http:https://julialang.org/license
+
 # these could fail on an embedded installation
 # but for now, we don't handle that case
 dlls = Libdl.dllist()

test/linalg/uniformscaling.jl

@@ -0,0 +1,92 @@
+# This file is a part of Julia. License is MIT: http:https://julialang.org/license
+
+using Base.Test
+
+srand(123)
+
+# Uniform scaling
+@test I[1,1] == 1 # getindex
+@test I[1,2] == 0 # getindex
+@test I === I' # transpose
+@test one(UniformScaling{Float32}) == UniformScaling(one(Float32))
+@test zero(UniformScaling{Float32}) == UniformScaling(zero(Float32))
+@test zero(UniformScaling(rand(Complex128))) == zero(UniformScaling{Complex128})
+@test one(UniformScaling(rand(Complex128))) == one(UniformScaling{Complex128})
+@test eltype(one(UniformScaling(rand(Complex128)))) == Complex128
+@test -one(UniformScaling(2)) == UniformScaling(-1)
+
+α = randn()
+@test α .* UniformScaling(1.0) == UniformScaling(1.0) .* α
+@test UniformScaling(α)./α == UniformScaling(1.0)
+@test α + UniformScaling(1.0) == UniformScaling(1.0) + α
+@test α - UniformScaling(1.0) == -(UniformScaling(1.0) - α)
+@test copy(UniformScaling(one(Float64))) == UniformScaling(one(Float64))
+@test sprint(show,UniformScaling(one(Float32))) == "Base.LinAlg.UniformScaling{Float32}\n1.0*I"
+
+λ = complex(randn(),randn())
+J = UniformScaling(λ)
+@test ndims(J) == 2
+@test transpose(J) == J
+@test J*eye(2) == conj(J'eye(2)) # ctranpose (and A(c)_mul_B)
+@test I + I === UniformScaling(2) # +
+@test inv(I) == I
+@test inv(J) == UniformScaling(inv(λ))
+
+B = bitrand(2,2)
+@test B + I == B + eye(B)
+@test I + B == B + eye(B)
+
+A = randn(2,2)
+@test A + I == A + eye(A)
+@test I + A == A + eye(A)
+@test I - I === UniformScaling(0)
+@test B - I == B - eye(B)
+@test I - B == eye(B) - B
+@test A - I == A - eye(A)
+@test I - A == eye(A) - A
+@test I*J === UniformScaling(λ)
+@test B*J == B*λ
+@test J*B == B*λ
+
+S = sprandn(3,3,0.5)
+@test S*J == S*λ
+@test J*S == S*λ
+@test A*J == A*λ
+@test J*A == A*λ
+@test J*ones(3) == ones(3)*λ
+@test λ*J === UniformScaling(λ*J.λ)
+@test J*λ === UniformScaling(λ*J.λ)
+@test J/I === J
+@test I/A == inv(A)
+@test A/I == A
+@test I/λ === UniformScaling(1/λ)
+@test I\J === J
+
+T = LowerTriangular(randn(3,3))
+@test T + J == full(T) + J
+@test J + T == J + full(T)
+@test T - J == full(T) - J
+@test J - T == J - full(T)
+@test T\I == inv(T)
+T = LinAlg.UnitLowerTriangular(randn(3,3))
+@test T + J == full(T) + J
+@test J + T == J + full(T)
+@test T - J == full(T) - J
+@test J - T == J - full(T)
+@test T\I == inv(T)
+T = UpperTriangular(randn(3,3))
+@test T + J == full(T) + J
+@test J + T == J + full(T)
+@test T - J == full(T) - J
+@test J - T == J - full(T)
+@test T\I == inv(T)
+T = LinAlg.UnitUpperTriangular(randn(3,3))
+@test T + J == full(T) + J
+@test J + T == J + full(T)
+@test T - J == full(T) - J
+@test J - T == J - full(T)
+@test T\I == inv(T)
+
+@test I\A == A
+@test A\I == inv(A)
+@test λ\I === UniformScaling(1/λ)

test/linalg2.jl

@@ -334,65 +334,6 @@ end
 @test_approx_eq norm(Float64[1e-300, 1], -3)*1e300 1
 @test_approx_eq norm(Float64[1e300, 1], 3)*1e-300 1
 
-# Uniform scaling
-@test I[1,1] == 1 # getindex
-@test I[1,2] == 0 # getindex
-@test I === I' # transpose
-@test one(UniformScaling{Float32}) == UniformScaling(one(Float32))
-@test zero(UniformScaling{Float32}) == UniformScaling(zero(Float32))
-@test zero(UniformScaling(rand(Complex128))) == zero(UniformScaling{Complex128})
-@test one(UniformScaling(rand(Complex128))) == one(UniformScaling{Complex128})
-@test eltype(one(UniformScaling(rand(Complex128)))) == Complex128
-@test -one(UniformScaling(2)) == UniformScaling(-1)
-α = randn()
-@test α .* UniformScaling(1.0) == UniformScaling(1.0) .* α
-@test UniformScaling(α)./α == UniformScaling(1.0)
-@test α + UniformScaling(1.0) == UniformScaling(1.0) + α
-@test α - UniformScaling(1.0) == -(UniformScaling(1.0) - α)
-@test copy(UniformScaling(one(Float64))) == UniformScaling(one(Float64))
-@test sprint(show,UniformScaling(one(Float32))) == "Base.LinAlg.UniformScaling{Float32}\n1.0*I"
-
-λ = complex(randn(),randn())
-J = UniformScaling(λ)
-@test ndims(J) == 2
-@test transpose(J) == J
-@test J*eye(2) == conj(J'eye(2)) # ctranpose (and A(c)_mul_B)
-@test I + I === UniformScaling(2) # +
-@test inv(I) == I
-@test inv(J) == UniformScaling(inv(λ))
-B = bitrand(2,2)
-@test B + I == B + eye(B)
-@test I + B == B + eye(B)
-A = randn(2,2)
-@test A + I == A + eye(A)
-@test I + A == A + eye(A)
-@test I - I === UniformScaling(0)
-@test B - I == B - eye(B)
-@test I - B == eye(B) - B
-@test A - I == A - eye(A)
-@test I - A == eye(A) - A
-@test I*J === UniformScaling(λ)
-@test B*J == B*λ
-@test J*B == B*λ
-S = sprandn(3,3,0.5)
-@test S*J == S*λ
-@test J*S == S*λ
-@test A*J == A*λ
-@test J*A == A*λ
-@test J*ones(3) == ones(3)*λ
-@test λ*J === UniformScaling(λ*J.λ)
-@test J*λ === UniformScaling(λ*J.λ)
-@test J/I === J
-@test I/A == inv(A)
-@test A/I == A
-@test I/λ === UniformScaling(1/λ)
-@test I\J === J
-T = LowerTriangular(randn(3,3))
-@test T\I == inv(T)
-@test I\A == A
-@test A\I == inv(A)
-@test λ\I === UniformScaling(1/λ)
-
 ## Issue related tests
 # issue #1447
 let

test/testdefs.jl

@@ -3,7 +3,7 @@
 using Base.Test
 
 function runtests(name)
- @printf(" \033[1m*\033[0m \033[31m%-20s\033[0m", name)
+ @printf(" \033[1m*\033[0m \033[31m%-21s\033[0m", name)
  tt = @elapsed include("$name.jl")
  @printf(" in %6.2f seconds\n", tt)
  nothing