Merge branch 'master' into SEBAtest

Author: mkalia94

src/LinearAlgebra.jl

@@ -13,7 +13,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     copy, copyto!, copymutable, cos, cosh, cot, coth, csc, csch, eltype, exp, fill!, floor,
     getindex, hcat, getproperty, imag, inv, invpermuterows!, isapprox, isequal, isone, iszero,
     IndexStyle, kron, kron!, length, log, map, ndims, one, oneunit, parent, permutecols!,
-    permutedims, permuterows!, power_by_squaring, promote_rule, real, sec, sech, setindex!,
+    permutedims, permuterows!, power_by_squaring, promote_rule, real, isreal, sec, sech, setindex!,
     show, similar, sin, sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc,
     typed_hcat, vec, view, zero
 import Base: AbstractArray, AbstractMatrix, Array, Matrix

src/bidiag.jl

@@ -286,6 +286,7 @@ axes(M::Bidiagonal) = (ax = axes(M.dv, 1); (ax, ax))
 for func in (:conj, :copy, :real, :imag)
     @eval ($func)(M::Bidiagonal) = Bidiagonal(($func)(M.dv), ($func)(M.ev), M.uplo)
 end
+isreal(M::Bidiagonal) = isreal(M.dv) && isreal(M.ev)
 
 adjoint(B::Bidiagonal{<:Number}) = Bidiagonal(vec(adjoint(B.dv)), vec(adjoint(B.ev)), B.uplo == 'U' ? :L : :U)
 adjoint(B::Bidiagonal{<:Number, <:Base.ReshapedArray{<:Number,1,<:Adjoint}}) =

src/diagonal.jl

@@ -256,6 +256,8 @@ factorize(D::Diagonal) = D
 real(D::Diagonal) = Diagonal(real(D.diag))
 imag(D::Diagonal) = Diagonal(imag(D.diag))
 
+isreal(D::Diagonal) = isreal(D.diag)
+
 iszero(D::Diagonal) = all(iszero, D.diag)
 isone(D::Diagonal) = all(isone, D.diag)
 isdiag(D::Diagonal) = all(isdiag, D.diag)

src/generic.jl

@@ -1358,8 +1358,15 @@ ishermitian(x::Number) = (x == conj(x))
 _iszero(V) = iszero(V)
 # A Base.FastContiguousSubArray view of a StridedArray
 FastContiguousSubArrayStrided{T,N,P<:StridedArray,I<:Tuple{AbstractUnitRange, Vararg{Any}}} = Base.SubArray{T,N,P,I,true}
-# using mapreduce instead of all permits vectorization
-_iszero(V::FastContiguousSubArrayStrided) = mapreduce(iszero, &, V, init=true)
+# Reducing over the entire array instead of calling `all` within `iszero` permits vectorization
+# The loop is equivalent to a mapreduce, but is faster to compile
+function _iszero(V::FastContiguousSubArrayStrided)
+    ret = true
+    for i in eachindex(V)
+        ret &= iszero(@inbounds V[i])
+    end
+    ret
+end
 
 """
     istriu(A::AbstractMatrix, k::Integer = 0) -> Bool

src/symmetric.jl

@@ -441,6 +441,10 @@ issymmetric(A::Hermitian{<:Real}) = true
 issymmetric(A::Hermitian{<:Complex}) = isreal(A)
 issymmetric(A::Symmetric) = true
 
+# check if the symmetry is known from the type
+_issymmetric(::Union{SymSymTri, Hermitian{<:Real}}) = true
+_issymmetric(::Any) = false
+
 adjoint(A::Hermitian) = A
 transpose(A::Symmetric) = A
 adjoint(A::Symmetric{<:Real}) = A

src/triangular.jl

@@ -376,6 +376,11 @@ istril(A::Transpose, k::Integer=0) = istriu(A.parent, -k)
 istriu(A::Adjoint, k::Integer=0) = istril(A.parent, -k)
 istriu(A::Transpose, k::Integer=0) = istril(A.parent, -k)
 
+istril(U::UpperTriangular, k::Integer=0) = istril(parent(U), max(-1, k))
+istril(U::UnitUpperTriangular, k::Integer=0) = k < 0 ? false : istril(parent(U), k)
+istriu(U::LowerTriangular, k::Integer=0) = istriu(parent(U), min(1, k))
+istriu(U::UnitLowerTriangular, k::Integer=0) = k > 0 ? false : istriu(parent(U), k)
+
 function tril!(A::UpperTriangular{T}, k::Integer=0) where {T}
     if k < 0
         fill!(A.data, zero(T))

src/tridiag.jl

@@ -111,16 +111,14 @@ function (::Type{SymTri})(A::AbstractMatrix) where {SymTri <: SymTridiagonal}
     checksquare(A)
     du = diag(A, 1)
     d  = diag(A)
-    dl = diag(A, -1)
-    if _checksymmetric(d, du, dl)
-        SymTri(d, du)
-    else
+    if !(_issymmetric(A) || _checksymmetric(d, du, diag(A, -1)))
         throw(ArgumentError("matrix is not symmetric; cannot convert to SymTridiagonal"))
     end
+    return SymTri(d, du)
 end
 
 _checksymmetric(d, du, dl) = all(((x, y),) -> x == transpose(y), zip(du, dl)) && all(issymmetric, d)
-_checksymmetric(A::AbstractMatrix) = _checksymmetric(diagview(A), diagview(A, 1), diagview(A, -1))
+_checksymmetric(A::AbstractMatrix) = _issymmetric(A) || _checksymmetric(diagview(A), diagview(A, 1), diagview(A, -1))
 
 SymTridiagonal{T,V}(S::SymTridiagonal{T,V}) where {T,V<:AbstractVector{T}} = S
 SymTridiagonal{T,V}(S::SymTridiagonal) where {T,V<:AbstractVector{T}} =
@@ -173,6 +171,7 @@ _copyto_banded!(dest::SymTridiagonal, src::SymTridiagonal) =
 for func in (:conj, :copy, :real, :imag)
     @eval ($func)(M::SymTridiagonal) = SymTridiagonal(($func)(M.dv), ($func)(M.ev))
 end
+isreal(S::SymTridiagonal) = isreal(S.dv) && isreal(S.ev)
 
 transpose(S::SymTridiagonal) = S
 adjoint(S::SymTridiagonal{<:Number}) = SymTridiagonal(vec(adjoint(S.dv)), vec(adjoint(S.ev)))
@@ -667,6 +666,7 @@ for func in (:conj, :copy, :real, :imag)
         Tridiagonal(($func)(M.dl), ($func)(M.d), ($func)(M.du))
     end
 end
+isreal(T::Tridiagonal) = isreal(T.dl) && isreal(T.d) && isreal(T.du)
 
 adjoint(S::Tridiagonal{<:Number}) = Tridiagonal(vec(adjoint(S.du)), vec(adjoint(S.d)), vec(adjoint(S.dl)))
 adjoint(S::Tridiagonal{<:Number, <:Base.ReshapedArray{<:Number,1,<:Adjoint}}) =

test/bidiag.jl

@@ -1172,4 +1172,12 @@ end
     @test_throws InexactError convert(Bidiagonal, M)
 end
 
+@testset "isreal" begin
+    M = Bidiagonal(ones(2), ones(1), :U)
+    @test @inferred((M -> Val(isreal(M)))(M)) == Val(true)
+    M = complex.(M)
+    @test isreal(M)
+    @test !isreal(im*M)
+end
+
 end # module TestBidiagonal

test/diagonal.jl

@@ -1481,4 +1481,12 @@ end
     @test ishermitian(D) == ishermitian(A)
 end
 
+@testset "isreal" begin
+    D = Diagonal(ones(2))
+    @test @inferred((D -> Val(isreal(D)))(D)) == Val(true)
+    D = complex.(D)
+    @test isreal(D)
+    @test !isreal(im*D)
+end
+
 end # module TestDiagonal

test/generic.jl

@@ -167,6 +167,12 @@ end
                 @test_throws DimensionMismatch mul!(similar(a), a, Diagonal(Vector{Float64}(undef, an+1)))
                 @test mul!(similar(a), a, Diagonal(1.:an)) == a.*Vector(1:an)'
                 @test mul!(similar(a), a, Diagonal(1:an))  == a.*Vector(1:an)'
+
+                @testset "different axes" begin
+                    O = OffsetArray(similar(a), ntuple(_->2, ndims(a)))
+                    @test mul!(O, a, 2) == OffsetArray(2a, axes(O))
+                    @test mul!(O, 2, a) == OffsetArray(2a, axes(O))
+                end
             end
 
             @testset "Scaling with 5-argument mul!" begin
@@ -180,6 +186,20 @@ end
                 @test mul!(copy(a), Diagonal([1; 2]), a, 10, 100)   == 10a.*[1; 2] .+ 100a
                 @test mul!(copy(a), a, Diagonal(1.:an), 10, 100) == 10a.*Vector(1:an)' .+ 100a
                 @test mul!(copy(a), a, Diagonal(1:an), 10, 100)  == 10a.*Vector(1:an)' .+ 100a
+
+                @testset "different axes" begin
+                    if eltype(a) <: Number
+                        O = OffsetArray(ones(size(a)), ntuple(_->2, ndims(a)))
+                        @test mul!(copy(O), a, 2, 3, 4) == OffsetArray(6a .+ 4, axes(O))
+                        @test mul!(copy(O), 2, a, 3, 4) == OffsetArray(6a .+ 4, axes(O))
+                        @test mul!(copy(O), a, 2, 3, 0) == OffsetArray(6a, axes(O))
+                        @test mul!(copy(O), 2, a, 3, 0) == OffsetArray(6a, axes(O))
+                        @test mul!(copy(O), a, 2, 1, 4) == OffsetArray(2a .+ 4, axes(O))
+                        @test mul!(copy(O), 2, a, 1, 4) == OffsetArray(2a .+ 4, axes(O))
+                        @test mul!(copy(O), a, 2, 1, 0) == OffsetArray(2a, axes(O))
+                        @test mul!(copy(O), 2, a, 1, 0) == OffsetArray(2a, axes(O))
+                    end
+                end
             end
         end
     end