From 84faacfef790c15af6bf1f499d6578a45973c79e Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Sat, 17 Nov 2018 18:27:03 +0100 Subject: [PATCH 1/4] proper diagonal in copytri! (fix #30055) --- stdlib/LinearAlgebra/src/matmul.jl | 5 +++++ stdlib/LinearAlgebra/src/symmetric.jl | 4 ++-- stdlib/LinearAlgebra/test/matmul.jl | 7 +++++++ 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 43459b1f9e4cf..6325fa5d3ccc5 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -340,6 +340,11 @@ function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false) else throw(ArgumentError("uplo argument must be 'U' (upper) or 'L' (lower), got $uplo")) end + if conjugate + for i = 1:n + A[i,i] = adjoint(A[i,i]) + end + end A end diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 0111b3faf4f93..f20e78d76dbb1 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -236,14 +236,14 @@ similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = s function Matrix(A::Symmetric) B = copytri!(convert(Matrix, copy(A.data)), A.uplo) for i = 1:size(A, 1) - B[i,i] = symmetric(B[i,i], sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) + B[i,i] = symmetric(A[i,i], sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) end return B end function Matrix(A::Hermitian) B = copytri!(convert(Matrix, copy(A.data)), A.uplo, true) for i = 1:size(A, 1) - B[i,i] = hermitian(B[i,i], sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) + B[i,i] = hermitian(A[i,i], sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) end return B end diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 62c7d7cc8faf2..1f0607798ca9c 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -301,6 +301,13 @@ end @test_throws ArgumentError LinearAlgebra.copytri!(Matrix{Float64}(undef,10,10),'Z') +@testset "Issue 30055" begin + A = UpperTriangular([1+im 2+im 3+im; 4+im 5+im 6+im; 7+im 9+im 0])' + @test copy(A) == A # copy uses copytri! + A = LowerTriangular([1+im 2+im 3+im; 4+im 5+im 6+im; 7+im 9+im 0])' + @test copy(A) == A # copy uses copytri! +end + @testset "gemv! and gemm_wrapper for $elty" for elty in [Float32,Float64,ComplexF64,ComplexF32] A10x10, x10, x11 = Array{elty}.(undef, ((10,10), 10, 11)) @test_throws DimensionMismatch LinearAlgebra.gemv!(x10,'N',A10x10,x11) From 7073d71458e8d88c01908c951a17361c395e0985 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Mon, 19 Nov 2018 14:06:37 +0100 Subject: [PATCH 2/4] added sprandn methods with Type --- stdlib/SparseArrays/src/sparsematrix.jl | 4 +++- stdlib/SparseArrays/src/sparsevector.jl | 2 ++ stdlib/SparseArrays/test/sparse.jl | 7 +++++++ stdlib/SparseArrays/test/sparsevector.jl | 5 +++++ 4 files changed, 17 insertions(+), 1 deletion(-) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 02c975168c6ec..7fac079ac97d6 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1471,7 +1471,7 @@ sprand(r::AbstractRNG, ::Type{Bool}, m::Integer, n::Integer, density::AbstractFl sprand(::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where {T} = sprand(GLOBAL_RNG, T, m, n, density) """ - sprandn([rng], m[,n],p::AbstractFloat) + sprandn([rng],[,Type],m[,n],p::AbstractFloat) Create a random sparse vector of length `m` or sparse matrix of size `m` by `n` with the specified (independent) probability `p` of any entry being nonzero, @@ -1488,6 +1488,8 @@ julia> sprandn(2, 2, 0.75) """ sprandn(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat) = sprand(r,m,n,density,randn,Float64) sprandn(m::Integer, n::Integer, density::AbstractFloat) = sprandn(GLOBAL_RNG,m,n,density) +sprandn(r::AbstractRNG, ::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where T = sprand(r,m,n,density,(r,i) -> randn(r,T,i), T) +sprandn(::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where T = sprandn(GLOBAL_RNG,T,m,n,density) LinearAlgebra.fillstored!(S::SparseMatrixCSC, x) = (fill!(nzvalview(S), x); S) diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index e7f0e1a2db26a..ce9fbbc35e947 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -506,6 +506,8 @@ sprand(::Type{T}, n::Integer, p::AbstractFloat) where {T} = sprand(GLOBAL_RNG, T sprandn(n::Integer, p::AbstractFloat) = sprand(GLOBAL_RNG, n, p, randn) sprandn(r::AbstractRNG, n::Integer, p::AbstractFloat) = sprand(r, n, p, randn) +sprandn(::Type{T}, n::Integer, p::AbstractFloat) where T = sprand(GLOBAL_RNG, n, p, (r, i) -> randn(r, T, i)) +sprandn(r::AbstractRNG, ::Type{T}, n::Integer, p::AbstractFloat) where T = sprand(r, n, p, (r, i) -> randn(r, T, i)) ## Indexing into Matrices can return SparseVectors diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index c8e5ebd8b5acb..9eb2069d3984c 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -2335,4 +2335,11 @@ end @test m2.module == SparseArrays end +@testset "sprandn with type $T" for T in (Float64, Float32, Float16, ComplexF64, ComplexF32, ComplexF16) + @test sprandn(T, 5, 5, 0.5) isa AbstractSparseMatrix{T} +end +@testset "sprandn with invalid type $T" for T in (AbstractFloat, BigFloat, Complex) + @test_throws MethodError sprandn(T, 5, 5, 0.5) +end + end # module diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index 6061b21a6721c..d8ada3963e2c5 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -154,6 +154,11 @@ end end end + let xr = sprandn(ComplexF64, 1000, 0.9) + @test isa(xr, SparseVector{ComplexF64,Int}) + @test length(xr) == 1000 + end + let xr = sprand(Bool, 1000, 0.9) @test isa(xr, SparseVector{Bool,Int}) @test length(xr) == 1000 From 49f415640139c3c8d59cb86aa0bf46e1f89feee9 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Mon, 7 Jan 2019 19:55:09 +0100 Subject: [PATCH 3/4] additional parameter in copytri! for diagonal --- stdlib/LinearAlgebra/src/matmul.jl | 16 ++++++++-------- stdlib/LinearAlgebra/src/triangular.jl | 16 ++++++++-------- stdlib/LinearAlgebra/test/matmul.jl | 22 ++++++++++++++++++---- 3 files changed, 34 insertions(+), 20 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 6325fa5d3ccc5..a752a1a5b6931 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -327,24 +327,24 @@ function mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, transB end # Supporting functions for matrix multiplication -function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false) +# copy transposed(adjoint) of upper(lower) side-digonals. Optionally include diagonal. +function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false, diag::Bool=false) + copytri!(A, Val(uplo), Val(conjugate), Val(diag)) +end +function copytri!(A::AbstractMatrix, ::Val{uplo}, ::Val{conjugate}, ::Val{diag}) where {uplo, conjugate, diag} n = checksquare(A) + off = diag ? 0 : 1 if uplo == 'U' - for i = 1:(n-1), j = (i+1):n + for i = 1:n, j = (i+off):n A[j,i] = conjugate ? adjoint(A[i,j]) : transpose(A[i,j]) end elseif uplo == 'L' - for i = 1:(n-1), j = (i+1):n + for i = 1:n, j = (i+off):n A[i,j] = conjugate ? adjoint(A[j,i]) : transpose(A[j,i]) end else throw(ArgumentError("uplo argument must be 'U' (upper) or 'L' (lower), got $uplo")) end - if conjugate - for i = 1:n - A[i,i] = adjoint(A[i,i]) - end - end A end diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index 5627cbb332d5e..0c66885d816cb 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -346,14 +346,14 @@ Base.copy(A::Transpose{<:Any,<:UpperTriangular}) = transpose!(copy(A.parent)) Base.copy(A::Transpose{<:Any,<:UnitLowerTriangular}) = transpose!(copy(A.parent)) Base.copy(A::Transpose{<:Any,<:UnitUpperTriangular}) = transpose!(copy(A.parent)) -transpose!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L')) -transpose!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L')) -transpose!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U')) -transpose!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U')) -adjoint!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L' , true)) -adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , true)) -adjoint!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U' , true)) -adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true)) +transpose!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L', false, true)) +transpose!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L', false, true)) +transpose!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U', false, true)) +transpose!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U', false, true)) +adjoint!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L' , true, true)) +adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , true, true)) +adjoint!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U' , true, true)) +adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true, true)) diag(A::LowerTriangular) = diag(A.data) diag(A::UnitLowerTriangular) = fill(one(eltype(A)), size(A,1)) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 1f0607798ca9c..bf929c3ba3b12 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -302,10 +302,24 @@ end @test_throws ArgumentError LinearAlgebra.copytri!(Matrix{Float64}(undef,10,10),'Z') @testset "Issue 30055" begin - A = UpperTriangular([1+im 2+im 3+im; 4+im 5+im 6+im; 7+im 9+im 0])' - @test copy(A) == A # copy uses copytri! - A = LowerTriangular([1+im 2+im 3+im; 4+im 5+im 6+im; 7+im 9+im 0])' - @test copy(A) == A # copy uses copytri! + B = [1+im 2+im 3+im; 4+im 5+im 6+im; 7+im 9+im im] + A = UpperTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' + A = LowerTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' + B = Matrix{Matrix{Complex{Int}}}(undef, 2, 2) + B[1,1] = [1+im 2+im; 3+im 4+im] + B[2,1] = [1+2im 1+3im;1+3im 1+4im] + B[1,2] = [7+im 8+2im; 9+3im 4im] + B[2,2] = [9+im 8+im; 7+im 6+im] + A = UpperTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' + A = LowerTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' end @testset "gemv! and gemm_wrapper for $elty" for elty in [Float32,Float64,ComplexF64,ComplexF32] From fa3d7ad5cfbf7f260e1381b5d9964098e6c723cc Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Tue, 8 Jan 2019 13:02:34 +0100 Subject: [PATCH 4/4] @inline copytri! to enforce constant propagation --- stdlib/LinearAlgebra/src/matmul.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index a752a1a5b6931..bbdf30fedaaea 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -328,10 +328,7 @@ end # Supporting functions for matrix multiplication # copy transposed(adjoint) of upper(lower) side-digonals. Optionally include diagonal. -function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false, diag::Bool=false) - copytri!(A, Val(uplo), Val(conjugate), Val(diag)) -end -function copytri!(A::AbstractMatrix, ::Val{uplo}, ::Val{conjugate}, ::Val{diag}) where {uplo, conjugate, diag} +@inline function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false, diag::Bool=false) n = checksquare(A) off = diag ? 0 : 1 if uplo == 'U'