From e62193f5f1daae67184aa2979b4da81b7c7ac863 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 15 Feb 2025 16:06:43 +0530 Subject: [PATCH 01/28] Remove `LinearAlgebra` qualifications in `cholesky.jl` (cherry picked from commit 95d009b49d78c9638a22d486eb23c31351e3f386) --- src/cholesky.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cholesky.jl b/src/cholesky.jl index 03f7c273..2952503e 100644 --- a/src/cholesky.jl +++ b/src/cholesky.jl @@ -305,7 +305,7 @@ function _cholpivoted!(A::AbstractMatrix, ::Type{UpperTriangular}, tol::Real, ch rTA = real(eltype(A)) # checks Base.require_one_based_indexing(A) - n = LinearAlgebra.checksquare(A) + n = checksquare(A) # initialization piv = collect(1:n) dots = zeros(rTA, n) @@ -354,7 +354,7 @@ function _cholpivoted!(A::AbstractMatrix, ::Type{LowerTriangular}, tol::Real, ch rTA = real(eltype(A)) # checks Base.require_one_based_indexing(A) - n = LinearAlgebra.checksquare(A) + n = checksquare(A) # initialization piv = collect(1:n) dots = zeros(rTA, n) @@ -813,7 +813,7 @@ function rdiv!(B::AbstractMatrix, C::Cholesky) end end -function LinearAlgebra.rdiv!(B::AbstractMatrix, C::CholeskyPivoted) +function rdiv!(B::AbstractMatrix, C::CholeskyPivoted) n = size(C, 2) for i in 1:size(B, 1) permute!(view(B, i, 1:n), C.piv) @@ -965,7 +965,7 @@ function lowrankdowndate!(C::Cholesky, v::AbstractVector) s = conj(v[i]/Aii) s2 = abs2(s) if s2 > 1 - throw(LinearAlgebra.PosDefException(i)) + throw(PosDefException(i)) end c = sqrt(1 - abs2(s)) From 0fa47f673865d50119952fd06da9465847238164 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 27 Feb 2025 17:12:50 +0530 Subject: [PATCH 02/28] Avoid materializing `diag` in `Diagonal` `kron` (#1230) (cherry picked from commit 907a2020c851e0067cf52ea1dfd4b443eb4ead9f) --- src/diagonal.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index c4f88276..ddee8981 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -743,16 +743,16 @@ end kron(A::Diagonal, B::Diagonal) = Diagonal(kron(A.diag, B.diag)) function kron(A::Diagonal, B::SymTridiagonal) - kdv = kron(diag(A), B.dv) + kdv = kron(A.diag, B.dv) # We don't need to drop the last element - kev = kron(diag(A), _pushzero(_evview(B))) + kev = kron(A.diag, _pushzero(_evview(B))) SymTridiagonal(kdv, kev) end function kron(A::Diagonal, B::Tridiagonal) # `_droplast!` is only guaranteed to work with `Vector` - kd = convert(Vector, kron(diag(A), B.d)) - kdl = _droplast!(convert(Vector, kron(diag(A), _pushzero(B.dl)))) - kdu = _droplast!(convert(Vector, kron(diag(A), _pushzero(B.du)))) + kd = convert(Vector, kron(A.diag, B.d)) + kdl = _droplast!(convert(Vector, kron(A.diag, _pushzero(B.dl)))) + kdu = _droplast!(convert(Vector, kron(A.diag, _pushzero(B.du)))) Tridiagonal(kdl, kd, kdu) end From b8258de2783e129089e6cf3da17db9a66491893f Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 17 Mar 2025 00:56:53 +0530 Subject: [PATCH 03/28] Reduce `stable_muladdmul` branches in `generic matvecmul!` (#1240) Since `beta` is not used in this operation, we may set it to a constant to eliminate two of the branches generated by `@stable_muladdmul`. (cherry picked from commit 5122dbf42ee83da56429eabe55de5b2a73344a75) --- src/matmul.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/matmul.jl b/src/matmul.jl index 84a83d3a..6eb9e9df 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -962,7 +962,7 @@ function __generic_matvecmul!(::typeof(identity), C::AbstractVector, A::Abstract end for k = eachindex(B) aoffs = (k-1)*Astride - b = @stable_muladdmul MulAddMul(alpha,beta)(B[k]) + b = @stable_muladdmul MulAddMul(alpha,false)(B[k]) for i = eachindex(C) C[i] += A[aoffs + i] * b end From d455b8387f78cc547a9fa491fe645bbbb9461a45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateus=20Ara=C3=BAjo?= Date: Tue, 1 Apr 2025 15:00:39 +0200 Subject: [PATCH 04/28] fix dispatch to herk (#1247) (cherry picked from commit e64a3df85ea640ca769d0e0065e8e1648dfc2bfe) --- src/matmul.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/matmul.jl b/src/matmul.jl index 6eb9e9df..bf3699be 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -692,7 +692,7 @@ end # the aggressive constprop pushes tA and tB into gemm_wrapper!, which is needed for wrap calls within it # to be concretely inferred Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T}, - alpha::Number, beta::Number) where {T<:BlasFloat} + α::Number, β::Number) where {T<:BlasFloat} nC = checksquare(C) tA_uc = uppercase(tA) # potentially convert a WrapperChar to a Char if tA_uc == 'T' @@ -708,8 +708,8 @@ Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::Abst # BLAS.syrk! only updates symmetric C # alternatively, make non-zero β a show-stopper for BLAS.syrk! - if iszero(beta) || issymmetric(C) - α, β = promote(alpha, beta, zero(T)) + if iszero(β) || issymmetric(C) + alpha, beta = promote(α, β, zero(T)) if (alpha isa Union{Bool,T} && beta isa Union{Bool,T} && stride(A, 1) == stride(C, 1) == 1 && @@ -717,7 +717,7 @@ Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::Abst return copytri!(BLAS.syrk!('U', tA, alpha, A, beta, C), 'U') end end - return gemm_wrapper!(C, tA, tAt, A, A, alpha, beta) + return gemm_wrapper!(C, tA, tAt, A, A, α, β) end # legacy method syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T}, _add::MulAddMul = MulAddMul()) where {T<:BlasFloat} = @@ -743,7 +743,7 @@ Base.@constprop :aggressive function herk_wrapper!(C::Union{StridedMatrix{T}, St # Result array does not need to be initialized as long as beta==0 # C = Matrix{T}(undef, mA, mA) - if iszero(β) || issymmetric(C) + if iszero(β) || ishermitian(C) alpha, beta = promote(α, β, zero(T)) if (alpha isa Union{Bool,T} && beta isa Union{Bool,T} && From 5e2e77701d04cebfffa6ff5aaeba5929f4ff0e4c Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Tue, 1 Apr 2025 18:56:15 -0400 Subject: [PATCH 05/28] use smaller matrix size in `peakflops` on 32-bit (cherry picked from commit 0c87d0d195129b5c65a9a27a46bde08d3c803ea8) --- src/LinearAlgebra.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/LinearAlgebra.jl b/src/LinearAlgebra.jl index 2c76dcdc..1526df5e 100644 --- a/src/LinearAlgebra.jl +++ b/src/LinearAlgebra.jl @@ -726,6 +726,8 @@ end (\)(F::TransposeFactorization{T,<:LU}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = ldiv(F, B) +const default_peakflops_size = Int === Int32 ? 2048 : 4096 + """ LinearAlgebra.peakflops(n::Integer=4096; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false) @@ -750,7 +752,7 @@ of the problem that is solved on each processor. This function requires at least Julia 1.1. In Julia 1.0 it is available from the standard library `InteractiveUtils`. """ -function peakflops(n::Integer=4096; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false) +function peakflops(n::Integer=default_peakflops_size; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false) t = zeros(Float64, ntrials) for i=1:ntrials a = ones(eltype,n,n) From 1efe04addb95056d85c4e3ab96ad465b2b896175 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 26 Apr 2025 14:14:27 +0530 Subject: [PATCH 06/28] Only `@noinline` error path in `matmul_size_check` (#1310) Instead of `@noinline` on the entire size-check function, we now separate the error-throwing part into a separate function and mark it as `@noinline`. This way, the size check may still be evaluated inline, and only the error path will not be inlined. This improves performance for small matmul. ```julia julia> A = [1 2; 3 4]; julia> @btime $A * $A; 53.361 ns (2 allocations: 112 bytes) # v"1.13.0-DEV.438" 47.504 ns (2 allocations: 112 bytes) # this PR ``` (cherry picked from commit e4e8c1967c6fd2a5b7a991d43265c861b9a9cae7) --- src/matmul.jl | 76 ++++++++++++++++++++++++++++----------------------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/src/matmul.jl b/src/matmul.jl index bf3699be..10a4c81c 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -412,52 +412,60 @@ lmul!(A, B) _vec_or_mat_str(s::Tuple{Any}) = "vector" _vec_or_mat_str(s::Tuple{Any,Any}) = "matrix" -@noinline function matmul_size_check(sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}}) +function matmul_size_check(sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}}) szA2 = get(sizeA, 2, 1) if szA2 != sizeB[1] - strA = _vec_or_mat_str(sizeA) - strB = _vec_or_mat_str(sizeB) - B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB - size_or_len_str_B = B_size_len isa Integer ? "length" : "size" - dim_or_len_str_B = B_size_len isa Integer ? "length" : "first dimension" - pos_str_A = LazyString(length(sizeA) == length(sizeB) ? "first " : "", strA) - pos_str_B = LazyString(length(sizeA) == length(sizeB) ? "second " : "", strB) - throw(DimensionMismatch( - LazyString( - lazy"incompatible dimensions for matrix multiplication: ", - lazy"tried to multiply a $strA of size $sizeA with a $strB of $size_or_len_str_B $B_size_len. ", - lazy"The second dimension of the $pos_str_A: $szA2, does not match the $dim_or_len_str_B of the $pos_str_B: $(sizeB[1])." - ) - ) - ) + matmul_size_check_error(sizeA, sizeB) end return nothing end -@noinline function matmul_size_check(sizeC::Tuple{Integer,Vararg{Integer}}, sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}}) +@noinline function matmul_size_check_error(sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}}) + strA = _vec_or_mat_str(sizeA) + strB = _vec_or_mat_str(sizeB) + szA2 = get(sizeA, 2, 1) + B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB + size_or_len_str_B = B_size_len isa Integer ? "length" : "size" + dim_or_len_str_B = B_size_len isa Integer ? "length" : "first dimension" + pos_str_A = LazyString(length(sizeA) == length(sizeB) ? "first " : "", strA) + pos_str_B = LazyString(length(sizeA) == length(sizeB) ? "second " : "", strB) + throw(DimensionMismatch( + LazyString( + "incompatible dimensions for matrix multiplication: ", + lazy"tried to multiply a $strA of size $sizeA with a $strB of $size_or_len_str_B $B_size_len. ", + lazy"The second dimension of the $pos_str_A: $szA2, does not match the $dim_or_len_str_B of the $pos_str_B: $(sizeB[1])." + ) + ) + ) +end +function matmul_size_check(sizeC::Tuple{Integer,Vararg{Integer}}, sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}}) matmul_size_check(sizeA, sizeB) szB2 = get(sizeB, 2, 1) szC2 = get(sizeC, 2, 1) if sizeC[1] != sizeA[1] || szC2 != szB2 - strA = _vec_or_mat_str(sizeA) - strB = _vec_or_mat_str(sizeB) - strC = _vec_or_mat_str(sizeC) - C_size_len = length(sizeC) == 1 ? sizeC[1] : sizeC - size_or_len_str_C = C_size_len isa Integer ? "length" : "size" - B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB - size_or_len_str_B = B_size_len isa Integer ? "length" : "size" - destsize = length(sizeB) == length(sizeC) == 1 ? sizeA[1] : (sizeA[1], szB2) - size_or_len_str_dest = destsize isa Integer ? "length" : "size" - throw(DimensionMismatch( - LazyString( - "incompatible destination size: ", - lazy"the destination $strC of $size_or_len_str_C $C_size_len is incomatible with the product of a $strA of size $sizeA and a $strB of $size_or_len_str_B $B_size_len. ", - lazy"The destination must be of $size_or_len_str_dest $destsize." - ) - ) - ) + matmul_size_check_error(sizeC, sizeA, sizeB) end return nothing end +@noinline function matmul_size_check_error(sizeC::Tuple{Integer,Vararg{Integer}}, sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}}) + strA = _vec_or_mat_str(sizeA) + strB = _vec_or_mat_str(sizeB) + strC = _vec_or_mat_str(sizeC) + szB2 = get(sizeB, 2, 1) + C_size_len = length(sizeC) == 1 ? sizeC[1] : sizeC + size_or_len_str_C = C_size_len isa Integer ? "length" : "size" + B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB + size_or_len_str_B = B_size_len isa Integer ? "length" : "size" + destsize = length(sizeB) == length(sizeC) == 1 ? sizeA[1] : (sizeA[1], szB2) + size_or_len_str_dest = destsize isa Integer ? "length" : "size" + throw(DimensionMismatch( + LazyString( + "incompatible destination size: ", + lazy"the destination $strC of $size_or_len_str_C $C_size_len is incomatible with the product of a $strA of size $sizeA and a $strB of $size_or_len_str_B $B_size_len. ", + lazy"The destination must be of $size_or_len_str_dest $destsize." + ) + ) + ) +end # We may inline the matmul2x2! and matmul3x3! calls for `α == true` # to simplify the @stable_muladdmul branches From a7eea43908489a5a317050ad995692e566096a82 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 28 Apr 2025 10:35:22 +0530 Subject: [PATCH 07/28] Refine column ranges in `_isbanded_impl` (#1267) Fixes the column ranges, so that these correctly correspond to the ranges over which there are top and bottom rows to consider beyond the bands. The change is to the upper limit of the range. (cherry picked from commit 95703b590aa06a4ab833650a4e9d0b98956474c7) --- src/generic.jl | 22 +++++++++++++++++++++- test/generic.jl | 21 ++++++++++++++++----- 2 files changed, 37 insertions(+), 6 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index 2b03b249..fc0b4664 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1486,6 +1486,24 @@ function _isbanded_impl(A, kl, ku) beyond ku, where the elements should all be zero. The reason we separate this from the third group is that we may loop over all the rows using A[:, col] instead of A[rowrange, col], which is usually faster. + + E.g., in the following 6x10 matrix with (kl,ku) = (-1,1): + 1 1 0 0 0 0 0 0 0 0 + 1 2 2 0 0 0 0 0 0 0 + 0 2 3 3 0 0 0 0 0 0 + 0 0 3 4 4 0 0 0 0 0 + 0 0 0 4 5 5 0 0 0 0 + 0 0 0 0 5 6 6 0 0 0 + + last_col_nonzeroblocks: 7, as every column beyond this is entirely zero + last_col_emptytoprows: 2, as there are zeros above the stored bands beyond this column + last_col_nonemptybottomrows: 4, as there are no zeros below the stored bands beyond this column + colrange_onlybottomrows: 1:2, as these columns only have zeros below the stored bands + colrange_topbottomrows: 3:4, as these columns have zeros both above and below the stored bands + colrange_onlytoprows_nonzero: 5:7, as these columns only have zeros above the stored bands + colrange_zero_block: 8:10, as every column in this range is filled with zeros + + These are used to determine which rows to check for zeros in each column. =# last_col_nonzeroblocks = size(A,1) + ku # fully zero rectangular block beyond this column @@ -1493,7 +1511,9 @@ function _isbanded_impl(A, kl, ku) last_col_nonemptybottomrows = size(A,1) + kl - 1 # empty bottom rows after this column colrange_onlybottomrows = firstindex(A,2):min(last_col_nonemptybottomrows, last_col_emptytoprows) - colrange_topbottomrows = max(last_col_emptytoprows, last(colrange_onlybottomrows))+1:last_col_nonzeroblocks + col_topbotrows_start = max(last_col_emptytoprows, last(colrange_onlybottomrows))+1 + col_topbotrows_end = min(last_col_nonemptybottomrows, last_col_nonzeroblocks) + colrange_topbottomrows = col_topbotrows_start:col_topbotrows_end colrange_onlytoprows_nonzero = last(colrange_topbottomrows)+1:last_col_nonzeroblocks colrange_zero_block = last_col_nonzeroblocks+1:lastindex(A,2) diff --git a/test/generic.jl b/test/generic.jl index 6d11ec82..fcbe4234 100644 --- a/test/generic.jl +++ b/test/generic.jl @@ -514,17 +514,17 @@ end @testset "generic functions for checking whether matrices have banded structure" begin pentadiag = [1 2 3; 4 5 6; 7 8 9] - tridiag = [1 2 0; 4 5 6; 0 8 9] - tridiagG = GenericArray([1 2 0; 4 5 6; 0 8 9]) + tridiag = diagm(-1=>1:6, 1=>1:6) + tridiagG = GenericArray(tridiag) Tridiag = Tridiagonal(tridiag) ubidiag = [1 2 0; 0 5 6; 0 0 9] - ubidiagG = GenericArray([1 2 0; 0 5 6; 0 0 9]) + ubidiagG = GenericArray(ubidiag) uBidiag = Bidiagonal(ubidiag, :U) lbidiag = [1 0 0; 4 5 0; 0 8 9] - lbidiagG = GenericArray([1 0 0; 4 5 0; 0 8 9]) + lbidiagG = GenericArray(lbidiag) lBidiag = Bidiagonal(lbidiag, :L) adiag = [1 0 0; 0 5 0; 0 0 9] - adiagG = GenericArray([1 0 0; 0 5 0; 0 0 9]) + adiagG = GenericArray(adiag) aDiag = Diagonal(adiag) @testset "istriu" begin @test !istriu(pentadiag) @@ -618,6 +618,17 @@ end end end end + + tridiag = diagm(-1=>1:6, 1=>1:6) + A = [tridiag zeros(size(tridiag,1), 2)] + G = GenericArray(A) + @testset for (kl,ku) in Iterators.product(-10:10, -10:10) + @test isbanded(A, kl, ku) == isbanded(G, kl, ku) + end + @testset for k in -10:10 + @test istriu(A,k) == istriu(G,k) + @test istril(A,k) == istril(G,k) + end end @testset "missing values" begin From c816ae201d89cae3dc7047664b1df83ba90f9aec Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 30 Apr 2025 23:27:41 +0530 Subject: [PATCH 08/28] Fix empty `Tridiagonal` broadcast (#1324) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fixes ```julia julia> T = Tridiagonal(1:0, 1:0, 1:0) 0×0 Tridiagonal{Int64, UnitRange{Int64}} julia> T .+ T ERROR: ArgumentError: invalid GenericMemory size: the number of elements is either negative or too large for system address width ``` After this, ```julia julia> T .+ T 0×0 Tridiagonal{Int64, Vector{Int64}} ``` The changes are minor, but the largish diff is because of adding a loop to the tests, and additional indentation as a consequence. (cherry picked from commit dccd6f810e4b69da91a208a57b1f796b4db31166) --- src/structuredbroadcast.jl | 4 +- test/structuredbroadcast.jl | 296 ++++++++++++++++++------------------ 2 files changed, 152 insertions(+), 148 deletions(-) diff --git a/src/structuredbroadcast.jl b/src/structuredbroadcast.jl index 2c1e7bfa..05684131 100644 --- a/src/structuredbroadcast.jl +++ b/src/structuredbroadcast.jl @@ -84,9 +84,9 @@ function structured_broadcast_alloc(bc, ::Type{Bidiagonal}, ::Type{ElType}, n) w return Bidiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n1), uplo) end structured_broadcast_alloc(bc, ::Type{SymTridiagonal}, ::Type{ElType}, n) where {ElType} = - SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n-1)) + SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, max(0,n-1))) structured_broadcast_alloc(bc, ::Type{Tridiagonal}, ::Type{ElType}, n) where {ElType} = - Tridiagonal(Array{ElType}(undef, n-1),Array{ElType}(undef, n),Array{ElType}(undef, n-1)) + Tridiagonal(Array{ElType}(undef, max(0,n-1)),Array{ElType}(undef, n),Array{ElType}(undef, max(0,n-1))) structured_broadcast_alloc(bc, ::Type{LowerTriangular}, ::Type{ElType}, n) where {ElType} = LowerTriangular(Array{ElType}(undef, n, n)) structured_broadcast_alloc(bc, ::Type{UpperTriangular}, ::Type{ElType}, n) where {ElType} = diff --git a/test/structuredbroadcast.jl b/test/structuredbroadcast.jl index 5b5cc0cd..a360dfdb 100644 --- a/test/structuredbroadcast.jl +++ b/test/structuredbroadcast.jl @@ -8,160 +8,164 @@ isdefined(Main, :SizedArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), using .Main.SizedArrays @testset "broadcast[!] over combinations of scalars, structured matrices, and dense vectors/matrices" begin - N = 10 - s = rand() - fV = rand(N) - fA = rand(N, N) - Z = copy(fA) - D = Diagonal(rand(N)) - B = Bidiagonal(rand(N), rand(N - 1), :U) - T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1)) - S = SymTridiagonal(rand(N), rand(N - 1)) - U = UpperTriangular(rand(N,N)) - L = LowerTriangular(rand(N,N)) - M = Matrix(rand(N,N)) - structuredarrays = (D, B, T, U, L, M, S) - fstructuredarrays = map(Array, structuredarrays) - for (X, fX) in zip(structuredarrays, fstructuredarrays) - @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX)) - @test broadcast!(sin, Z, X) == broadcast(sin, fX) - @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) - @test broadcast!(cos, Z, X) == broadcast(cos, fX) - @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX)) - @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) - @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) - @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) - @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) - @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) - - @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 - @test X .* 2.0 isa typeof(X) - @test X .* (2.0,) isa typeof(X) - @test isequal(X .* Inf, fX .* Inf) - - two = 2 - @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two - @test X .^ 2 isa typeof(X) - @test X .^ (2,) isa typeof(X) - @test X .^ two isa typeof(X) - @test X .^ 0 == fX .^ 0 - @test X .^ -1 == fX .^ -1 - - for (Y, fY) in zip(structuredarrays, fstructuredarrays) - @test broadcast(+, X, Y) == broadcast(+, fX, fY) - @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) - @test broadcast(*, X, Y) == broadcast(*, fX, fY) - @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + @testset for N in (0,1,2,10) # some edge cases, and a structured case + s = rand() + fV = rand(N) + fA = rand(N, N) + Z = copy(fA) + D = Diagonal(rand(N)) + B = Bidiagonal(rand(N), rand(max(0,N-1)), :U) + T = Tridiagonal(rand(max(0,N-1)), rand(N), rand(max(0,N-1))) + S = SymTridiagonal(rand(N), rand(max(0,N-1))) + U = UpperTriangular(rand(N,N)) + L = LowerTriangular(rand(N,N)) + M = Matrix(rand(N,N)) + structuredarrays = (D, B, T, U, L, M, S) + fstructuredarrays = map(Array, structuredarrays) + @testset "$(nameof(typeof(X)))" for (X, fX) in zip(structuredarrays, fstructuredarrays) + @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX)) + @test broadcast!(sin, Z, X) == broadcast(sin, fX) + @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) + @test broadcast!(cos, Z, X) == broadcast(cos, fX) + @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX)) + @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) + @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) + @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) + @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) + @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) + + @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 + @test X .* 2.0 isa typeof(X) + @test X .* (2.0,) isa typeof(X) + @test isequal(X .* Inf, fX .* Inf) + + two = 2 + @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two + @test X .^ 2 isa typeof(X) + @test X .^ (2,) isa typeof(X) + @test X .^ two isa typeof(X) + @test X .^ 0 == fX .^ 0 + @test X .^ -1 == fX .^ -1 + + for (Y, fY) in zip(structuredarrays, fstructuredarrays) + @test broadcast(+, X, Y) == broadcast(+, fX, fY) + @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) + @test broadcast(*, X, Y) == broadcast(*, fX, fY) + @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + end end - end - diagonals = (D, B, T) - fdiagonals = map(Array, diagonals) - for (X, fX) in zip(diagonals, fdiagonals) - for (Y, fY) in zip(diagonals, fdiagonals) - @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY) - @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) - @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY) - @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + diagonals = (D, B, T) + fdiagonals = map(Array, diagonals) + for (X, fX) in zip(diagonals, fdiagonals) + for (Y, fY) in zip(diagonals, fdiagonals) + @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY) + @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) + @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY) + @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + end end - end - UU = UnitUpperTriangular(rand(N,N)) - UL = UnitLowerTriangular(rand(N,N)) - unittriangulars = (UU, UL) - Ttris = typeof.((UpperTriangular(parent(UU)), LowerTriangular(parent(UU)))) - funittriangulars = map(Array, unittriangulars) - for (X, fX, Ttri) in zip(unittriangulars, funittriangulars, Ttris) - @test (Q = broadcast(sin, X); typeof(Q) == Ttri && Q == broadcast(sin, fX)) - @test broadcast!(sin, Z, X) == broadcast(sin, fX) - @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) - @test broadcast!(cos, Z, X) == broadcast(cos, fX) - @test (Q = broadcast(*, s, X); typeof(Q) == Ttri && Q == broadcast(*, s, fX)) - @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) - @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) - @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) - @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) - @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) - - @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 - @test X .* 2.0 isa Ttri - @test X .* (2.0,) isa Ttri - @test isequal(X .* Inf, fX .* Inf) - - two = 2 - @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two - @test X .^ 2 isa typeof(X) # special cased, as isstructurepreserving - @test X .^ (2,) isa Ttri - @test X .^ two isa Ttri - @test X .^ 0 == fX .^ 0 - @test X .^ -1 == fX .^ -1 - - for (Y, fY) in zip(unittriangulars, funittriangulars) - @test broadcast(+, X, Y) == broadcast(+, fX, fY) - @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) - @test broadcast(*, X, Y) == broadcast(*, fX, fY) - @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + UU = UnitUpperTriangular(rand(N,N)) + UL = UnitLowerTriangular(rand(N,N)) + unittriangulars = (UU, UL) + Ttris = typeof.((UpperTriangular(parent(UU)), LowerTriangular(parent(UU)))) + funittriangulars = map(Array, unittriangulars) + for (X, fX, Ttri) in zip(unittriangulars, funittriangulars, Ttris) + @test (Q = broadcast(sin, X); typeof(Q) == Ttri && Q == broadcast(sin, fX)) + @test broadcast!(sin, Z, X) == broadcast(sin, fX) + @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) + @test broadcast!(cos, Z, X) == broadcast(cos, fX) + @test (Q = broadcast(*, s, X); typeof(Q) == Ttri && Q == broadcast(*, s, fX)) + @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) + @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) + @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) + @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) + @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) + + @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 + @test X .* 2.0 isa Ttri + @test X .* (2.0,) isa Ttri + @test isequal(X .* Inf, fX .* Inf) + + two = 2 + @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two + @test X .^ 2 isa typeof(X) # special cased, as isstructurepreserving + @test X .^ (2,) isa Ttri + @test X .^ two isa Ttri + @test X .^ 0 == fX .^ 0 + @test X .^ -1 == fX .^ -1 + + for (Y, fY) in zip(unittriangulars, funittriangulars) + @test broadcast(+, X, Y) == broadcast(+, fX, fY) + @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) + @test broadcast(*, X, Y) == broadcast(*, fX, fY) + @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + end end - end - @testset "type-stability in Bidiagonal" begin - B2 = @inferred (B -> .- B)(B) - @test B2 isa Bidiagonal - @test B2 == -1 * B - B2 = @inferred (B -> B .* 2)(B) - @test B2 isa Bidiagonal - @test B2 == B + B - B2 = @inferred (B -> 2 .* B)(B) - @test B2 isa Bidiagonal - @test B2 == B + B - B2 = @inferred (B -> B ./ 1)(B) - @test B2 isa Bidiagonal - @test B2 == B - B2 = @inferred (B -> 1 .\ B)(B) - @test B2 isa Bidiagonal - @test B2 == B + @testset "type-stability in Bidiagonal" begin + B2 = @inferred (B -> .- B)(B) + @test B2 isa Bidiagonal + @test B2 == -1 * B + B2 = @inferred (B -> B .* 2)(B) + @test B2 isa Bidiagonal + @test B2 == B + B + B2 = @inferred (B -> 2 .* B)(B) + @test B2 isa Bidiagonal + @test B2 == B + B + B2 = @inferred (B -> B ./ 1)(B) + @test B2 isa Bidiagonal + @test B2 == B + B2 = @inferred (B -> 1 .\ B)(B) + @test B2 isa Bidiagonal + @test B2 == B + end end end @testset "broadcast! where the destination is a structured matrix" begin - N = 5 - A = rand(N, N) - sA = A + copy(A') - D = Diagonal(rand(N)) - Bu = Bidiagonal(rand(N), rand(N - 1), :U) - Bl = Bidiagonal(rand(N), rand(N - 1), :L) - T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1)) - ◣ = LowerTriangular(rand(N,N)) - ◥ = UpperTriangular(rand(N,N)) - M = Matrix(rand(N,N)) - - @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D)) - @test broadcast!(sin, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U) - @test broadcast!(sin, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L) - @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T)) - @test broadcast!(sin, copy(◣), ◣) == LowerTriangular(sin.(◣)) - @test broadcast!(sin, copy(◥), ◥) == UpperTriangular(sin.(◥)) - @test broadcast!(sin, copy(M), M) == Matrix(sin.(M)) - @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A)) - @test broadcast!(*, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U) - @test broadcast!(*, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L) - @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A)) - @test broadcast!(*, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A)) - @test broadcast!(*, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A)) - @test broadcast!(*, copy(M), M, A) == Matrix(broadcast(*, M, A)) - - @test_throws ArgumentError broadcast!(cos, copy(D), D) == Diagonal(sin.(D)) - @test_throws ArgumentError broadcast!(cos, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U) - @test_throws ArgumentError broadcast!(cos, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L) - @test_throws ArgumentError broadcast!(cos, copy(T), T) == Tridiagonal(sin.(T)) - @test_throws ArgumentError broadcast!(cos, copy(◣), ◣) == LowerTriangular(sin.(◣)) - @test_throws ArgumentError broadcast!(cos, copy(◥), ◥) == UpperTriangular(sin.(◥)) - @test_throws ArgumentError broadcast!(+, copy(D), D, A) == Diagonal(broadcast(*, D, A)) - @test_throws ArgumentError broadcast!(+, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U) - @test_throws ArgumentError broadcast!(+, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L) - @test_throws ArgumentError broadcast!(+, copy(T), T, A) == Tridiagonal(broadcast(*, T, A)) - @test_throws ArgumentError broadcast!(+, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A)) - @test_throws ArgumentError broadcast!(+, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A)) - @test_throws ArgumentError broadcast!(*, copy(◥), ◣, 2) - @test_throws ArgumentError broadcast!(*, copy(Bu), Bl, 2) + @testset for N in (0,1,2,5) + A = rand(N, N) + sA = A + copy(A') + D = Diagonal(rand(N)) + Bu = Bidiagonal(rand(N), rand(max(0,N-1)), :U) + Bl = Bidiagonal(rand(N), rand(max(0,N-1)), :L) + T = Tridiagonal(rand(max(0,N-1)), rand(N), rand(max(0,N-1))) + ◣ = LowerTriangular(rand(N,N)) + ◥ = UpperTriangular(rand(N,N)) + M = Matrix(rand(N,N)) + + @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D)) + @test broadcast!(sin, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U) + @test broadcast!(sin, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L) + @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T)) + @test broadcast!(sin, copy(◣), ◣) == LowerTriangular(sin.(◣)) + @test broadcast!(sin, copy(◥), ◥) == UpperTriangular(sin.(◥)) + @test broadcast!(sin, copy(M), M) == Matrix(sin.(M)) + @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A)) + @test broadcast!(*, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U) + @test broadcast!(*, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L) + @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A)) + @test broadcast!(*, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A)) + @test broadcast!(*, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A)) + @test broadcast!(*, copy(M), M, A) == Matrix(broadcast(*, M, A)) + + if N > 2 + @test_throws ArgumentError broadcast!(cos, copy(D), D) + @test_throws ArgumentError broadcast!(cos, copy(Bu), Bu) + @test_throws ArgumentError broadcast!(cos, copy(Bl), Bl) + @test_throws ArgumentError broadcast!(cos, copy(T), T) + @test_throws ArgumentError broadcast!(cos, copy(◣), ◣) + @test_throws ArgumentError broadcast!(cos, copy(◥), ◥) + @test_throws ArgumentError broadcast!(+, copy(D), D, A) + @test_throws ArgumentError broadcast!(+, copy(Bu), Bu, A) + @test_throws ArgumentError broadcast!(+, copy(Bl), Bl, A) + @test_throws ArgumentError broadcast!(+, copy(T), T, A) + @test_throws ArgumentError broadcast!(+, copy(◣), ◣, A) + @test_throws ArgumentError broadcast!(+, copy(◥), ◥, A) + @test_throws ArgumentError broadcast!(*, copy(◥), ◣, 2) + @test_throws ArgumentError broadcast!(*, copy(Bu), Bl, 2) + end + end end @testset "map[!] over combinations of structured matrices" begin From 60f547e96e7d2ca91343a4a6a5355f8ebed0482f Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 1 May 2025 09:46:02 +0530 Subject: [PATCH 09/28] `iszero` check in hessenberg setindex (#1327) Instead of explicitly checking that the value `== 0`, we should ideally call `iszero` on it. (cherry picked from commit 9d139f5a15f505f9cb816b6f04a1a30c1294c0a8) --- src/hessenberg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 056b97e3..158e81c8 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -103,8 +103,8 @@ Base._reverse(A::UpperHessenberg, dims) = reverse!(Matrix(A); dims) Base.@propagate_inbounds function setindex!(A::UpperHessenberg, x, i::Integer, j::Integer) if i > j+1 - x == 0 || throw(ArgumentError("cannot set index in the lower triangular part " * - lazy"($i, $j) of an UpperHessenberg matrix to a nonzero value ($x)")) + iszero(x) || throw(ArgumentError(LazyString("cannot set index in the lower triangular part ", + lazy"($i, $j) of an UpperHessenberg matrix to a nonzero value ($x)"))) else A.data[i,j] = x end From 6d36322b948b75885ae886e97c484ddbd88d7023 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Sun, 4 May 2025 08:04:33 -0400 Subject: [PATCH 10/28] add missing methods for division of Hessenberg matrices (#1322) As noted [on discourse](https://discourse.julialang.org/t/add-a-upperhessenberg-row-in-linearalgebra-factorize/128530/4?u=stevengj), we were missing `/` and `\` methods for `UpperHessenberg` matrices, despite the existence of optimized `ldiv!` and `rdiv!` methods, so it was falling back to slow $O(n^3)$ methods. While I was at it, I also added methods for transpose/adjoint Hessenberg matrices. --------- Co-authored-by: Jishnu Bhattacharya (cherry picked from commit 9d26faf6b7346bbfb8ab935920dcaeee2367b7f8) --- src/hessenberg.jl | 23 +++++++++++++++++++++++ test/hessenberg.jl | 7 +++++++ 2 files changed, 30 insertions(+) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 158e81c8..297ea737 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -177,6 +177,29 @@ function \(U::UnitUpperTriangular, H::UpperHessenberg) UpperHessenberg(HH) end +function (\)(H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}}, B::AbstractVecOrMat) + TFB = typeof(oneunit(eltype(H)) \ oneunit(eltype(B))) + return ldiv!(H, copy_similar(B, TFB)) +end + +function (/)(B::AbstractMatrix, H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}}) + TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(H))) + return rdiv!(copy_similar(B, TFB), H) +end + +ldiv!(H::AdjOrTrans{<:Any,<:UpperHessenberg}, B::AbstractVecOrMat) = + (rdiv!(wrapperop(H)(B), parent(H)); B) +rdiv!(B::AbstractVecOrMat, H::AdjOrTrans{<:Any,<:UpperHessenberg}) = + (ldiv!(parent(H), wrapperop(H)(B)); B) + +# fix method ambiguities for right division, from adjtrans.jl: +/(u::AdjointAbsVec, A::UpperHessenberg) = adjoint(adjoint(A) \ u.parent) +/(u::TransposeAbsVec, A::UpperHessenberg) = transpose(transpose(A) \ u.parent) +/(u::AdjointAbsVec, A::Adjoint{<:Any,<:UpperHessenberg}) = adjoint(adjoint(A) \ u.parent) +/(u::TransposeAbsVec, A::Transpose{<:Any,<:UpperHessenberg}) = transpose(transpose(A) \ u.parent) +/(u::AdjointAbsVec, A::Transpose{<:Any,<:UpperHessenberg}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent) +/(u::TransposeAbsVec, A::Adjoint{<:Any,<:UpperHessenberg}) = transpose(conj(A.parent) \ u.parent) + # Solving (H+µI)x = b: we can do this in O(m²) time and O(m) memory # (in-place in x) by the RQ algorithm from: # diff --git a/test/hessenberg.jl b/test/hessenberg.jl index ebab3f29..f12ff7df 100644 --- a/test/hessenberg.jl +++ b/test/hessenberg.jl @@ -64,6 +64,13 @@ let n = 10 H = UpperHessenberg(Areal) @test Array(Hc + H) == Array(Hc) + Array(H) @test Array(Hc - H) == Array(Hc) - Array(H) + @testset "ldiv and rdiv" begin + for b in (b_, B_), H in (H, Hc, H', Hc', transpose(Hc)) + @test H * (H \ b) ≈ b + @test (b' / H) * H ≈ (Matrix(b') / H) * H ≈ b' + @test (transpose(b) / H) * H ≈ (Matrix(transpose(b)) / H) * H ≈ transpose(b) + end + end @testset "Preserve UpperHessenberg shape (issue #39388)" begin H = UpperHessenberg(Areal) A = rand(n,n) From 1ee6ec1e04b42171425a54d9ed1122a5c25a214d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 4 May 2025 21:25:38 +0530 Subject: [PATCH 11/28] Unwrap triangular matrices in broadcast (#1332) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In broadcasting over triangular matrices, we loop only over the stored indices. For these indices, indexing into a triangular matrix is equivalent to indexing into its parent. We may therefore replace an `UpperOrLowerTriangular` matrix by its parent, which removes the branch in `getindex`. This improves performance: ```julia julia> L = LowerTriangular(zeros(600,600)); julia> L2 = copy(L); julia> @btime broadcast!(+, $L2, $L, $L); 161.176 μs (0 allocations: 0 bytes) # master 80.894 μs (0 allocations: 0 bytes) # this PR ``` This replacement is performed recursively on a `Broadcasted` object by looping over its `args`, and non-triangular elements are left untouched. Only `UpperOrLowerTriangular` matrices will be replaced by their `parent`s. (cherry picked from commit 5165fd3f733c2651271de0c55defae9fc8ec7ce3) --- src/structuredbroadcast.jl | 16 ++++++++++++++-- src/triangular.jl | 6 ++---- test/structuredbroadcast.jl | 8 ++++++++ 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/src/structuredbroadcast.jl b/src/structuredbroadcast.jl index 05684131..cef4a788 100644 --- a/src/structuredbroadcast.jl +++ b/src/structuredbroadcast.jl @@ -264,13 +264,24 @@ function copyto!(dest::Tridiagonal, bc::Broadcasted{<:StructuredMatrixStyle}) return dest end +# Recursively replace wrapped matrices by their parents to improve broadcasting performance +# We may do this because the indexing within `copyto!` is restricted to the stored indices +preprocess_broadcasted(::Type{T}, A) where {T} = _preprocess_broadcasted(T, A) +function preprocess_broadcasted(::Type{T}, bc::Broadcasted) where {T} + args = map(x -> preprocess_broadcasted(T, x), bc.args) + Broadcast.Broadcasted(bc.f, args, bc.axes) +end +_preprocess_broadcasted(::Type{LowerTriangular}, A) = lowertridata(A) +_preprocess_broadcasted(::Type{UpperTriangular}, A) = uppertridata(A) + function copyto!(dest::LowerTriangular, bc::Broadcasted{<:StructuredMatrixStyle}) isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc)) axs = axes(dest) axes(bc) == axs || Broadcast.throwdm(axes(bc), axs) + bc_unwrapped = preprocess_broadcasted(LowerTriangular, bc) for j in axs[2] for i in j:axs[1][end] - @inbounds dest.data[i,j] = bc[CartesianIndex(i, j)] + @inbounds dest.data[i,j] = bc_unwrapped[CartesianIndex(i, j)] end end return dest @@ -280,9 +291,10 @@ function copyto!(dest::UpperTriangular, bc::Broadcasted{<:StructuredMatrixStyle} isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc)) axs = axes(dest) axes(bc) == axs || Broadcast.throwdm(axes(bc), axs) + bc_unwrapped = preprocess_broadcasted(UpperTriangular, bc) for j in axs[2] for i in 1:j - @inbounds dest.data[i,j] = bc[CartesianIndex(i, j)] + @inbounds dest.data[i,j] = bc_unwrapped[CartesianIndex(i, j)] end end return dest diff --git a/src/triangular.jl b/src/triangular.jl index 4a03cff6..1e750467 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -658,10 +658,8 @@ end uppertridata(A) = A lowertridata(A) = A -# we restrict these specializations only to strided matrices to avoid cases where an UpperTriangular type -# doesn't share its indexing with the parent -uppertridata(A::UpperTriangular{<:Any, <:StridedMatrix}) = parent(A) -lowertridata(A::LowerTriangular{<:Any, <:StridedMatrix}) = parent(A) +uppertridata(A::UpperTriangular) = parent(A) +lowertridata(A::LowerTriangular) = parent(A) @inline _rscale_add!(A::AbstractTriangular, B::AbstractTriangular, C::Number, alpha::Number, beta::Number) = @stable_muladdmul _triscale!(A, B, C, MulAddMul(alpha, beta)) diff --git a/test/structuredbroadcast.jl b/test/structuredbroadcast.jl index a360dfdb..641c761d 100644 --- a/test/structuredbroadcast.jl +++ b/test/structuredbroadcast.jl @@ -385,4 +385,12 @@ end @test ind == CartesianIndex(1,1) end +@testset "nested triangular broadcast" begin + for T in (LowerTriangular, UpperTriangular) + L = T(rand(Int,4,4)) + M = Matrix(L) + @test L .+ L .+ 0 .+ L .+ 0 .- L == 2M + end +end + end From f4bf551a503ffe78b4d101d5046903aa600a45cd Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 5 May 2025 22:24:50 +0530 Subject: [PATCH 12/28] Change `1:size` to `axes` in bidiag mul (#1337) This makes snippets of the code easier to follow (cherry picked from commit 2702a49494cfa2ba38c40098c47cbbb68f3b4e74) --- src/bidiag.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/bidiag.jl b/src/bidiag.jl index 89e5711c..ce5b2d80 100644 --- a/src/bidiag.jl +++ b/src/bidiag.jl @@ -963,7 +963,7 @@ function _mul_bitrisym!(C::AbstractVecOrMat, A::Bidiagonal, B::AbstractVecOrMat, if A.uplo == 'U' u = A.ev @inbounds begin - for j = 1:nB + for j in axes(B,2) b₀, b₊ = B[1, j], B[2, j] _modify!(_add, d[1]*b₀ + u[1]*b₊, C, (1, j)) for i = 2:nA - 1 @@ -976,7 +976,7 @@ function _mul_bitrisym!(C::AbstractVecOrMat, A::Bidiagonal, B::AbstractVecOrMat, else l = A.ev @inbounds begin - for j = 1:nB + for j in axes(B,2) b₀, b₊ = B[1, j], B[2, j] _modify!(_add, d[1]*b₀, C, (1, j)) for i = 2:nA - 1 @@ -996,7 +996,7 @@ function _mul_bitrisym!(C::AbstractVecOrMat, A::TriSym, B::AbstractVecOrMat, _ad d = _diag(A, 0) u = _diag(A, 1) @inbounds begin - for j = 1:nB + for j in axes(B,2) b₀, b₊ = B[1, j], B[2, j] _modify!(_add, d[1]*b₀ + u[1]*b₊, C, (1, j)) for i = 2:nA - 1 @@ -1028,7 +1028,7 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::TriSym, _add::MulAddMul) B21 = Bl[1] Bmm = Bd[m] Bm₋1m = Bu[m-1] - for i in 1:n + for i in axes(A,1) _modify!(_add, A[i,1] * B11 + A[i, 2] * B21, C, (i, 1)) _modify!(_add, A[i, m-1] * Bm₋1m + A[i, m] * Bmm, C, (i, m)) end @@ -1037,7 +1037,7 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::TriSym, _add::MulAddMul) Bj₋1j = Bu[j-1] Bjj = Bd[j] Bj₊1j = Bl[j] - for i = 1:n + for i in axes(A,1) _modify!(_add, A[i, j-1] * Bj₋1j + A[i, j]*Bjj + A[i, j+1] * Bj₊1j, C, (i, j)) end end @@ -1052,17 +1052,17 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::Bidiagonal, _add::MulAdd (iszero(m) || iszero(n)) && return C iszero(_add.alpha) && return _rmul_or_fill!(C, _add.beta) @inbounds if B.uplo == 'U' - for j in n:-1:2, i in 1:m + for j in reverse(axes(A,2)[2:end]), i in axes(A,1) _modify!(_add, A[i,j] * B.dv[j] + A[i,j-1] * B.ev[j-1], C, (i, j)) end - for i in 1:m + for i in axes(A,1) _modify!(_add, A[i,1] * B.dv[1], C, (i, 1)) end else # uplo == 'L' - for j in 1:n-1, i in 1:m + for j in axes(A,2)[1:end-1], i in axes(A,1) _modify!(_add, A[i,j] * B.dv[j] + A[i,j+1] * B.ev[j], C, (i, j)) end - for i in 1:m + for i in axes(A,1) _modify!(_add, A[i,n] * B.dv[n], C, (i, n)) end end @@ -1138,7 +1138,7 @@ function _dibimul!(C::AbstractMatrix, A::Diagonal, B::Bidiagonal, _add) if B.uplo == 'L' C[2,1] += _add(Ad[2]*Bev[1]) end - for col in 2:n-1 + for col in axes(A,1)[2:end-1] evrow = col+rowshift C[evrow, col] += _add(Ad[evrow]*Bev[col - evshift]) C[col, col] += _add(Ad[col]*Bdv[col]) @@ -1219,7 +1219,7 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector) @inbounds if B.uplo == 'U' x₀ = x[1] r = dot(x[1], dv[1], y[1]) - for j in 2:nx-1 + for j in eachindex(x)[2:end-1] x₋, x₀ = x₀, x[j] r += dot(adjoint(ev[j-1])*x₋ + adjoint(dv[j])*x₀, y[j]) end @@ -1229,7 +1229,7 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector) x₀ = x[1] x₊ = x[2] r = dot(adjoint(dv[1])*x₀ + adjoint(ev[1])*x₊, y[1]) - for j in 2:nx-1 + for j in eachindex(x)[2:end-1] x₀, x₊ = x₊, x[j+1] r += dot(adjoint(dv[j])*x₀ + adjoint(ev[j])*x₊, y[j]) end @@ -1260,15 +1260,15 @@ function ldiv!(c::AbstractVecOrMat, A::Bidiagonal, b::AbstractVecOrMat) zi = findfirst(iszero, A.dv) isnothing(zi) || throw(SingularException(zi)) - @inbounds for j in 1:nb + @inbounds for j in axes(b,2) if A.uplo == 'L' #do colwise forward substitution c[1,j] = bi1 = A.dv[1] \ b[1,j] - for i in 2:N + for i in eachindex(A.dv)[2:end] c[i,j] = bi1 = A.dv[i] \ (b[i,j] - A.ev[i - 1] * bi1) end else #do colwise backward substitution c[N,j] = bi1 = A.dv[N] \ b[N,j] - for i in (N - 1):-1:1 + for i in reverse(eachindex(A.ev)) c[i,j] = bi1 = A.dv[i] \ (b[i,j] - A.ev[i] * bi1) end end From 5dc7640ba1906b58cc7d3a8a47a4bd5d80405229 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 9 May 2025 10:49:10 +0530 Subject: [PATCH 13/28] `Char` uplo in `Bidiagonal` constructor (#1342) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This is an internal constructor, but many methods involving `Bidiagonal` matrices currently use the `Char` `uplo` directly to construct a new `Bidiagonal` matrix. We therefore may require the constructors to accept a `Char` as well, wherever they accept a `Symbol`. This already works for some of the `Bidiagonal` constructors, and this PR plugs a missing case. Fixes issues like ```julia julia> B = Bidiagonal(fill(SMatrix{2,3}(1:6), 4), fill(SMatrix{1,3}(1:3), 3), :L) 4×4 Bidiagonal{StaticArraysCore.SArray{S, Int64, 2} where S<:Tuple, Vector{StaticArraysCore.SArray{S, Int64, 2} where S<:Tuple}}: [1 3 5; 2 4 6] ⋅ ⋅ ⋅ [1 2 3] [1 3 5; 2 4 6] ⋅ ⋅ ⋅ [1 2 3] [1 3 5; 2 4 6] ⋅ ⋅ ⋅ [1 2 3] [1 3 5; 2 4 6] julia> 2B ERROR: MethodError: no method matching Bidiagonal(::Vector{SMatrix{2, 3, Int64, 6}}, ::Vector{SMatrix{1, 3, Int64, 3}}, ::Char) The type `Bidiagonal` exists, but no method is defined for this combination of argument types when trying to construct it. Closest candidates are: Bidiagonal(::Vector{T}, ::Vector{S}, ::Symbol) where {T, S} @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra_2.jl/src/bidiag.jl:77 Bidiagonal(::V, ::V, ::AbstractChar) where {T, V<:AbstractVector{T}} @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra_2.jl/src/bidiag.jl:71 Bidiagonal(::V, ::V, ::Symbol) where {T, V<:AbstractVector{T}} @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra_2.jl/src/bidiag.jl:68 ... ``` After this, ```julia julia> 2B 4×4 Bidiagonal{StaticArraysCore.SArray{S, Int64, 2} where S<:Tuple, Vector{StaticArraysCore.SArray{S, Int64, 2} where S<:Tuple}}: [2 6 10; 4 8 12] ⋅ ⋅ ⋅ [2 4 6] [2 6 10; 4 8 12] ⋅ ⋅ ⋅ [2 4 6] [2 6 10; 4 8 12] ⋅ ⋅ ⋅ [2 4 6] [2 6 10; 4 8 12] ``` We may require that each method passes a `Symbol` as `uplo`, but this would be a much larger change, and the current approach seems to be less disruptive. (cherry picked from commit 560276ba63d17654c7c709ae5e2ed04ba4dc74e1) --- src/bidiag.jl | 7 ++----- test/bidiag.jl | 4 ++++ 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/bidiag.jl b/src/bidiag.jl index ce5b2d80..ca95e386 100644 --- a/src/bidiag.jl +++ b/src/bidiag.jl @@ -65,16 +65,13 @@ julia> Bl = Bidiagonal(dv, ev, :L) # ev is on the first subdiagonal ⋅ ⋅ 9 4 ``` """ -function Bidiagonal(dv::V, ev::V, uplo::Symbol) where {T,V<:AbstractVector{T}} - Bidiagonal{T,V}(dv, ev, uplo) -end -function Bidiagonal(dv::V, ev::V, uplo::AbstractChar) where {T,V<:AbstractVector{T}} +function Bidiagonal(dv::V, ev::V, uplo::Union{Symbol, AbstractChar}) where {T,V<:AbstractVector{T}} Bidiagonal{T,V}(dv, ev, uplo) end #To allow Bidiagonal's where the "dv" is Vector{T} and "ev" Vector{S}, #where T and S can be promoted -function Bidiagonal(dv::Vector{T}, ev::Vector{S}, uplo::Symbol) where {T,S} +function Bidiagonal(dv::Vector{T}, ev::Vector{S}, uplo::Union{Symbol, AbstractChar}) where {T,S} TS = promote_type(T,S) return Bidiagonal{TS,Vector{TS}}(dv, ev, uplo) end diff --git a/test/bidiag.jl b/test/bidiag.jl index a39fa027..8d6a7d05 100644 --- a/test/bidiag.jl +++ b/test/bidiag.jl @@ -506,6 +506,10 @@ Random.seed!(1) @test Matrix{ComplexF64}(BD) == BD end +@testset "Constructors with Char uplo" begin + @test Bidiagonal(Int8[1,2], [1], 'U') == Bidiagonal(Int8[1,2], [1], :U) +end + # Issue 10742 and similar let A = Bidiagonal([1,2,3], [0,0], :U) @test istril(A) From aba980cea673156d35487f83da2bd04f64a29036 Mon Sep 17 00:00:00 2001 From: WalterMadelim <143272683+WalterMadelim@users.noreply.github.com> Date: Sun, 11 May 2025 20:58:26 +0800 Subject: [PATCH 14/28] Update the docstring of ldiv! (#1344) ![image](https://github.com/user-attachments/assets/724d0ca3-fc05-46a7-b22a-cc0ab929bfe0) A comment about the existing docstring is: The lower example reads counter-intuitive, _given_ the experience from reading the upper example. And `X` and `Y` are not like `A` and `B`. Hope the revised version reads clearer and intelligible. (cherry picked from commit 2d35f07f6c27d29cacc0c3f774032a6dbd96a0cb) --- src/LinearAlgebra.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/LinearAlgebra.jl b/src/LinearAlgebra.jl index 1526df5e..d2f82705 100644 --- a/src/LinearAlgebra.jl +++ b/src/LinearAlgebra.jl @@ -393,13 +393,13 @@ control over the factorization of `A`. ```jldoctest julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2]; -julia> X = [1; 2.5; 3]; +julia> B = [1, 2.5, 3]; -julia> Y = zero(X); +julia> Y = similar(B); # use similar since there is no need to read from it -julia> ldiv!(Y, qr(A), X); +julia> ldiv!(Y, qr(A), B); # you may also try qr!(A) to further reduce allocation -julia> Y ≈ A\\X +julia> Y ≈ A \\ B true ``` """ @@ -425,13 +425,13 @@ control over the factorization of `A`. ```jldoctest julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2]; -julia> X = [1; 2.5; 3]; +julia> B = [1, 2.5, 3]; -julia> Y = copy(X); +julia> B0 = copy(B); # a backup copy to facilitate testing -julia> ldiv!(qr(A), X); +julia> ldiv!(lu(A), B); # you may also try lu!(A) to further reduce allocation -julia> X ≈ A\\Y +julia> B ≈ A \\ B0 true ``` """ From 150e2fd4247f82806b22a473bd5d01357c6fed3e Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 11 May 2025 19:05:54 +0530 Subject: [PATCH 15/28] Fix scaling block unit triangular matrices (cherry picked from commit 2148fa45a68a4767a6f1663d7fd50fd2be70ef11) --- src/triangular.jl | 31 +++++++++++++------------------ test/triangular.jl | 10 ++++++++++ 2 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/triangular.jl b/src/triangular.jl index 1e750467..c4da1adb 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -1270,6 +1270,15 @@ end # Generic routines # #################### +function _set_diag!(B::UpperOrLowerTriangular, x) + # get a mutable array to modify the diagonal + Bm = parent(B) isa StridedArray ? B : copy!(similar(B), B) + for i in diagind(Bm.data, IndexStyle(Bm.data)) + Bm.data[i] = x + end + Bm +end + for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), (LowerTriangular, UnitLowerTriangular)) tstrided = t{<:Any, <:StridedMaybeAdjOrTransMat} @@ -1283,10 +1292,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), function (*)(A::$unitt, x::Number) B = $t(A.data)*x - for i in axes(A, 1) - B.data[i,i] = x - end - return B + _set_diag!(B, oneunit(eltype(A)) * x) end (*)(x::Number, A::$t) = $t(x*A.data) @@ -1298,10 +1304,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), function (*)(x::Number, A::$unitt) B = x*$t(A.data) - for i in axes(A, 1) - B.data[i,i] = x - end - return B + _set_diag!(B, x * oneunit(eltype(A))) end (/)(A::$t, x::Number) = $t(A.data/x) @@ -1313,11 +1316,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), function (/)(A::$unitt, x::Number) B = $t(A.data)/x - invx = inv(x) - for i in axes(A, 1) - B.data[i,i] = invx - end - return B + _set_diag!(B, oneunit(eltype(A)) / x) end (\)(x::Number, A::$t) = $t(x\A.data) @@ -1329,11 +1328,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), function (\)(x::Number, A::$unitt) B = x\$t(A.data) - invx = inv(x) - for i in axes(A, 1) - B.data[i,i] = invx - end - return B + _set_diag!(B, x \ oneunit(eltype(A))) end end end diff --git a/test/triangular.jl b/test/triangular.jl index eb6814b9..d63f8b06 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -923,4 +923,14 @@ end end end +@testset "block unit triangular scaling" begin + m = SizedArrays.SizedArray{(2,2)}([1 2; 3 4]) + U = UnitUpperTriangular(fill(m, 4, 4)) + M = Matrix{eltype(U)}(U) + @test U/2 == M/2 + @test 2\U == 2\M + @test U*2 == M*2 + @test 2*U == 2*M +end + end # module TestTriangular From 21e43466012225da8ae8b6e7f72ff6288114897e Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 12 May 2025 22:43:48 +0530 Subject: [PATCH 16/28] Revert "add missing methods for division of Hessenberg matrices (#1322)" This reverts commit 6d36322b948b75885ae886e97c484ddbd88d7023. --- src/hessenberg.jl | 23 ----------------------- test/hessenberg.jl | 7 ------- 2 files changed, 30 deletions(-) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 297ea737..158e81c8 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -177,29 +177,6 @@ function \(U::UnitUpperTriangular, H::UpperHessenberg) UpperHessenberg(HH) end -function (\)(H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}}, B::AbstractVecOrMat) - TFB = typeof(oneunit(eltype(H)) \ oneunit(eltype(B))) - return ldiv!(H, copy_similar(B, TFB)) -end - -function (/)(B::AbstractMatrix, H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}}) - TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(H))) - return rdiv!(copy_similar(B, TFB), H) -end - -ldiv!(H::AdjOrTrans{<:Any,<:UpperHessenberg}, B::AbstractVecOrMat) = - (rdiv!(wrapperop(H)(B), parent(H)); B) -rdiv!(B::AbstractVecOrMat, H::AdjOrTrans{<:Any,<:UpperHessenberg}) = - (ldiv!(parent(H), wrapperop(H)(B)); B) - -# fix method ambiguities for right division, from adjtrans.jl: -/(u::AdjointAbsVec, A::UpperHessenberg) = adjoint(adjoint(A) \ u.parent) -/(u::TransposeAbsVec, A::UpperHessenberg) = transpose(transpose(A) \ u.parent) -/(u::AdjointAbsVec, A::Adjoint{<:Any,<:UpperHessenberg}) = adjoint(adjoint(A) \ u.parent) -/(u::TransposeAbsVec, A::Transpose{<:Any,<:UpperHessenberg}) = transpose(transpose(A) \ u.parent) -/(u::AdjointAbsVec, A::Transpose{<:Any,<:UpperHessenberg}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent) -/(u::TransposeAbsVec, A::Adjoint{<:Any,<:UpperHessenberg}) = transpose(conj(A.parent) \ u.parent) - # Solving (H+µI)x = b: we can do this in O(m²) time and O(m) memory # (in-place in x) by the RQ algorithm from: # diff --git a/test/hessenberg.jl b/test/hessenberg.jl index f12ff7df..ebab3f29 100644 --- a/test/hessenberg.jl +++ b/test/hessenberg.jl @@ -64,13 +64,6 @@ let n = 10 H = UpperHessenberg(Areal) @test Array(Hc + H) == Array(Hc) + Array(H) @test Array(Hc - H) == Array(Hc) - Array(H) - @testset "ldiv and rdiv" begin - for b in (b_, B_), H in (H, Hc, H', Hc', transpose(Hc)) - @test H * (H \ b) ≈ b - @test (b' / H) * H ≈ (Matrix(b') / H) * H ≈ b' - @test (transpose(b) / H) * H ≈ (Matrix(transpose(b)) / H) * H ≈ transpose(b) - end - end @testset "Preserve UpperHessenberg shape (issue #39388)" begin H = UpperHessenberg(Areal) A = rand(n,n) From fdfa1a881aac7866994f1c8170e424d6936fe576 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 May 2025 12:19:21 +0530 Subject: [PATCH 17/28] Use method deletion instead of custom sysimage (#1312) (#1356) This backports #1312 to the backports-release-1.12 branch. Co-authored-by: Kristoffer Carlsson --- .buildkite/runtests.yml | 12 +- .ci/Manifest.toml | 222 --------------------------- .ci/Project.toml | 8 - .ci/create_sysimage.jl | 16 -- .ci/create_sysimage_and_run_docs.jl | 6 - .ci/create_sysimage_and_run_tests.jl | 7 - .ci/run_tests.jl | 11 ++ .github/workflows/ci.yml | 3 +- docs/make.jl | 2 + test/abstractq.jl | 2 + test/addmul.jl | 2 + test/adjtrans.jl | 2 + test/bidiag.jl | 2 + test/bitarray.jl | 2 + test/blas.jl | 2 + test/bunchkaufman.jl | 2 + test/cholesky.jl | 2 + test/dense.jl | 2 + test/diagonal.jl | 2 + test/eigen.jl | 2 + test/factorization.jl | 3 + test/generic.jl | 2 + test/givens.jl | 2 + test/hessenberg.jl | 4 +- test/lapack.jl | 2 + test/ldlt.jl | 2 + test/lq.jl | 2 + test/lu.jl | 2 + test/matmul.jl | 2 + test/pinv.jl | 2 + test/prune_old_LA.jl | 69 +++++++++ test/qr.jl | 2 + test/runtests.jl | 3 + test/schur.jl | 2 + test/special.jl | 2 + test/structuredbroadcast.jl | 3 + test/svd.jl | 2 + test/symmetric.jl | 2 + test/symmetriceigen.jl | 2 + test/triangular.jl | 2 + test/triangular2.jl | 2 + test/triangular3.jl | 2 + test/tridiag.jl | 2 + test/uniformscaling.jl | 2 + test/unitful.jl | 2 + 45 files changed, 161 insertions(+), 270 deletions(-) delete mode 100644 .ci/Manifest.toml delete mode 100644 .ci/Project.toml delete mode 100644 .ci/create_sysimage.jl delete mode 100644 .ci/create_sysimage_and_run_docs.jl delete mode 100644 .ci/create_sysimage_and_run_tests.jl create mode 100644 .ci/run_tests.jl create mode 100644 test/prune_old_LA.jl diff --git a/.buildkite/runtests.yml b/.buildkite/runtests.yml index a21eb633..936dc49d 100644 --- a/.buildkite/runtests.yml +++ b/.buildkite/runtests.yml @@ -16,8 +16,7 @@ steps: version: "nightly" arch: "i686" command: | - julia --color=yes --project=.ci -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project=.ci .ci/create_sysimage_and_run_tests.jl + julia --color=yes --code-coverage=@ .ci/run_tests.jl agents: queue: "julia" os: "linux" @@ -40,8 +39,7 @@ steps: - JuliaCI/julia#v1: version: "nightly" command: | - julia --color=yes --project=.ci --code-coverage=@ -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project=.ci --code-coverage=@ .ci/create_sysimage_and_run_tests.jl + julia --color=yes --code-coverage=@ .ci/run_tests.jl agents: queue: "julia" os: "linux" @@ -55,8 +53,7 @@ steps: - JuliaCI/julia-coverage#v1: codecov: true command: | - julia --color=yes --project=.ci --code-coverage=@ -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project=.ci --code-coverage=@ .ci/create_sysimage_and_run_tests.jl + julia --color=yes --code-coverage=@ .ci/run_tests.jl agents: queue: "julia" os: "macos" @@ -69,8 +66,7 @@ steps: - JuliaCI/julia-coverage#v1: codecov: true command: | - julia --color=yes --project=.ci --code-coverage=@ -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project=.ci --code-coverage=@ .ci/create_sysimage_and_run_tests.jl + julia --color=yes --code-coverage=@ .ci/run_tests.jl agents: queue: "julia" os: "windows" diff --git a/.ci/Manifest.toml b/.ci/Manifest.toml deleted file mode 100644 index 4f888b42..00000000 --- a/.ci/Manifest.toml +++ /dev/null @@ -1,222 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.12.0-DEV" -manifest_format = "2.0" -project_hash = "b25b51369ac1c4a63619c6f77655347661ec5439" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.2" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -version = "1.11.0" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" -version = "1.11.0" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.3.0+1" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" -version = "1.11.0" - -[[deps.Distributed]] -deps = ["Random", "Serialization", "Sockets"] -uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" -version = "1.11.0" - -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" -version = "1.11.0" - -[[deps.Glob]] -git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496" -uuid = "c27321d9-0574-5035-807b-f59d2c89b15c" -version = "1.3.1" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -version = "1.11.0" - -[[deps.JuliaSyntaxHighlighting]] -deps = ["StyledStrings"] -uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011" -version = "1.12.0" - -[[deps.LazyArtifacts]] -deps = ["Artifacts", "Pkg"] -uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" -version = "1.11.0" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.4" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.11.1+1" - -[[deps.LibGit2]] -deps = ["LibGit2_jll", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" -version = "1.11.0" - -[[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll"] -uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.9.0+0" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "OpenSSL_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.3+1" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -version = "1.11.0" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -path = ".." -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -version = "1.12.0" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" -version = "1.11.0" - -[[deps.Markdown]] -deps = ["Base64", "JuliaSyntaxHighlighting", "StyledStrings"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" -version = "1.11.0" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2024.12.31" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.3.0" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.29+0" - -[[deps.OpenSSL_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.15+2" - -[[deps.PackageCompiler]] -deps = ["Artifacts", "Glob", "LazyArtifacts", "Libdl", "Pkg", "Printf", "RelocatableFolders", "TOML", "UUIDs", "p7zip_jll"] -git-tree-sha1 = "5d13e5b70011762b74f86fc08385303589f80272" -uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" -version = "2.2.0" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.12.0" - - [deps.Pkg.extensions] - REPLExt = "REPL" - - [deps.Pkg.weakdeps] - REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" -version = "1.11.0" - -[[deps.Random]] -deps = ["SHA"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -version = "1.11.0" - -[[deps.RelocatableFolders]] -deps = ["SHA", "Scratch"] -git-tree-sha1 = "ffdaf70d81cf6ff22c2b6e733c900c3321cab864" -uuid = "05181044-ff0b-4ac5-8273-598c1e38db00" -version = "1.0.1" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -version = "0.7.0" - -[[deps.Scratch]] -deps = ["Dates"] -git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386" -uuid = "6c6a2e73-6563-6170-7368-637461726353" -version = "1.2.1" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" -version = "1.11.0" - -[[deps.Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" -version = "1.11.0" - -[[deps.StyledStrings]] -uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b" -version = "1.11.0" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - -[[deps.Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -version = "1.11.0" - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" -version = "1.11.0" - -[[deps.Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" -version = "1.11.0" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.3.1+2" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.12.0+0" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.64.0+1" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.5.0+2" diff --git a/.ci/Project.toml b/.ci/Project.toml deleted file mode 100644 index 3916fe72..00000000 --- a/.ci/Project.toml +++ /dev/null @@ -1,8 +0,0 @@ -[deps] -Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" -Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/.ci/create_sysimage.jl b/.ci/create_sysimage.jl deleted file mode 100644 index 43aea023..00000000 --- a/.ci/create_sysimage.jl +++ /dev/null @@ -1,16 +0,0 @@ -using PackageCompiler, Libdl - -function create_sysimage_linear_algebra() - sysimage = tempname() * "." * Libdl.dlext - - if haskey(ENV, "BUILDKITE") - ncores = Sys.CPU_THREADS - else - ncores = ceil(Int, Sys.CPU_THREADS / 2) - end - - withenv("JULIA_IMAGE_THREADS" => ncores) do - create_sysimage(["LinearAlgebra", "Test", "Distributed", "Dates", "Printf", "Random"]; sysimage_path=sysimage, incremental=false, filter_stdlibs=true) - end - return sysimage, ncores -end diff --git a/.ci/create_sysimage_and_run_docs.jl b/.ci/create_sysimage_and_run_docs.jl deleted file mode 100644 index 376219bb..00000000 --- a/.ci/create_sysimage_and_run_docs.jl +++ /dev/null @@ -1,6 +0,0 @@ -include("create_sysimage.jl") -sysimage, ncores = create_sysimage_linear_algebra() -doc_file = joinpath(@__DIR__, "..", "docs", "make.jl") -withenv("JULIA_NUM_THREADS" => 1) do - run(`$(Base.julia_cmd()) --sysimage=$sysimage --project=$(dirname(doc_file)) $doc_file`) -end diff --git a/.ci/create_sysimage_and_run_tests.jl b/.ci/create_sysimage_and_run_tests.jl deleted file mode 100644 index 071fe1fd..00000000 --- a/.ci/create_sysimage_and_run_tests.jl +++ /dev/null @@ -1,7 +0,0 @@ -include("create_sysimage.jl") -sysimage, ncores = create_sysimage_linear_algebra() -current_dir = @__DIR__ -cmd = """Base.runtests(["LinearAlgebra"]; propagate_project=true, ncores=$ncores)""" -withenv("JULIA_NUM_THREADS" => 1) do - run(`$(Base.julia_cmd()) --sysimage=$sysimage --project=$current_dir -e $cmd`) -end diff --git a/.ci/run_tests.jl b/.ci/run_tests.jl new file mode 100644 index 00000000..8a54ef62 --- /dev/null +++ b/.ci/run_tests.jl @@ -0,0 +1,11 @@ +if haskey(ENV, "BUILDKITE") + ncores = Sys.CPU_THREADS +else + ncores = ceil(Int, Sys.CPU_THREADS / 2) +end + +proj = abspath(joinpath(@__DIR__, "..")) +cmd = """Base.runtests(["LinearAlgebra"]; propagate_project=true, ncores=$ncores)""" +withenv("JULIA_NUM_THREADS" => 1) do + run(`$(Base.julia_cmd()) --project=$proj --compiled-modules=existing -e $cmd`) +end diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7c423f04..c1611fbf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,9 +21,8 @@ jobs: version: 'nightly' - name: Generate docs run: | - julia --color=yes --project=.ci -e 'using Pkg; Pkg.instantiate()' julia --color=yes --project=docs -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project=.ci .ci/create_sysimage_and_run_docs.jl + julia --color=yes --project=docs --compiled-modules=existing docs/make.jl env: DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/docs/make.jl b/docs/make.jl index 351fec43..a1d3a2d0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,3 +1,5 @@ +include("../test/prune_old_LA.jl") + using LinearAlgebra using Documenter: DocMeta, makedocs, deploydocs, HTML diff --git a/test/abstractq.jl b/test/abstractq.jl index 302f2abf..ff1499ae 100644 --- a/test/abstractq.jl +++ b/test/abstractq.jl @@ -2,6 +2,8 @@ module TestAbstractQ +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test using LinearAlgebra using LinearAlgebra: AbstractQ, AdjointQ diff --git a/test/addmul.jl b/test/addmul.jl index 903e3b17..f6b49d4e 100644 --- a/test/addmul.jl +++ b/test/addmul.jl @@ -2,6 +2,8 @@ module TestAddmul +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Base: rtoldefault using Test using LinearAlgebra diff --git a/test/adjtrans.jl b/test/adjtrans.jl index 50565583..bdeaa697 100644 --- a/test/adjtrans.jl +++ b/test/adjtrans.jl @@ -2,6 +2,8 @@ module TestAdjointTranspose +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") diff --git a/test/bidiag.jl b/test/bidiag.jl index 8d6a7d05..982eaf2d 100644 --- a/test/bidiag.jl +++ b/test/bidiag.jl @@ -2,6 +2,8 @@ module TestBidiagonal +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasReal, BlasFloat diff --git a/test/bitarray.jl b/test/bitarray.jl index 76fc21a6..a185873e 100644 --- a/test/bitarray.jl +++ b/test/bitarray.jl @@ -1,3 +1,5 @@ +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using LinearAlgebra, Test, Random tc(r1::NTuple{N,Any}, r2::NTuple{N,Any}) where {N} = all(x->tc(x...), [zip(r1,r2)...]) diff --git a/test/blas.jl b/test/blas.jl index ff1d8bab..b7f4a03a 100644 --- a/test/blas.jl +++ b/test/blas.jl @@ -2,6 +2,8 @@ module TestBLAS +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasReal, BlasComplex using Libdl: dlsym, dlopen diff --git a/test/bunchkaufman.jl b/test/bunchkaufman.jl index 68c519d1..cf4fa488 100644 --- a/test/bunchkaufman.jl +++ b/test/bunchkaufman.jl @@ -2,6 +2,8 @@ module TestBunchKaufman +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted using Base: getproperty diff --git a/test/cholesky.jl b/test/cholesky.jl index 6ba72432..4c1cb14b 100644 --- a/test/cholesky.jl +++ b/test/cholesky.jl @@ -2,6 +2,8 @@ module TestCholesky +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException, RankDeficientException, chkfullrank diff --git a/test/dense.jl b/test/dense.jl index 1d4b330c..47efb78c 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -2,6 +2,8 @@ module TestDense +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal using Test: GenericArray diff --git a/test/diagonal.jl b/test/diagonal.jl index d1f46f47..f83ad392 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -2,6 +2,8 @@ module TestDiagonal +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasFloat, BlasComplex diff --git a/test/eigen.jl b/test/eigen.jl index a82c7454..67e58d34 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -2,6 +2,8 @@ module TestEigen +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, UtiAUi! diff --git a/test/factorization.jl b/test/factorization.jl index f80c5197..b215c246 100644 --- a/test/factorization.jl +++ b/test/factorization.jl @@ -1,6 +1,9 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license module TestFactorization + +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra @testset "equality for factorizations - $f" for f in Any[ diff --git a/test/generic.jl b/test/generic.jl index fcbe4234..11ff874f 100644 --- a/test/generic.jl +++ b/test/generic.jl @@ -2,6 +2,8 @@ module TestGeneric +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using Test: GenericArray using LinearAlgebra: isbanded diff --git a/test/givens.jl b/test/givens.jl index 3726afa7..fb09ec4e 100644 --- a/test/givens.jl +++ b/test/givens.jl @@ -2,6 +2,8 @@ module TestGivens +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: Givens, Rotation, givensAlgorithm diff --git a/test/hessenberg.jl b/test/hessenberg.jl index ebab3f29..91f6d155 100644 --- a/test/hessenberg.jl +++ b/test/hessenberg.jl @@ -2,6 +2,8 @@ module TestHessenberg +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") @@ -283,5 +285,5 @@ end @test Q2'Q2 ≈ I end end - + end # module TestHessenberg diff --git a/test/lapack.jl b/test/lapack.jl index b4e3afe3..6e7bab0b 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -2,6 +2,8 @@ module TestLAPACK +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasInt diff --git a/test/ldlt.jl b/test/ldlt.jl index 51abf310..3454ef43 100644 --- a/test/ldlt.jl +++ b/test/ldlt.jl @@ -2,6 +2,8 @@ module TestLDLT +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random Random.seed!(123) diff --git a/test/lq.jl b/test/lq.jl index af7ab600..0d55f84f 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -2,6 +2,8 @@ module TestLQ +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul! diff --git a/test/lu.jl b/test/lu.jl index f5b4e247..ad4faa16 100644 --- a/test/lu.jl +++ b/test/lu.jl @@ -2,6 +2,8 @@ module TestLU +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: ldiv!, BlasReal, BlasInt, BlasFloat, rdiv! diff --git a/test/matmul.jl b/test/matmul.jl index 86c75ae5..d405225b 100644 --- a/test/matmul.jl +++ b/test/matmul.jl @@ -2,6 +2,8 @@ module TestMatmul +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Base: rtoldefault using Test, LinearAlgebra, Random using LinearAlgebra: mul!, Symmetric, Hermitian diff --git a/test/pinv.jl b/test/pinv.jl index c7268865..ecee397a 100644 --- a/test/pinv.jl +++ b/test/pinv.jl @@ -2,6 +2,8 @@ module TestPinv +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random Random.seed!(12345) diff --git a/test/prune_old_LA.jl b/test/prune_old_LA.jl new file mode 100644 index 00000000..f628cddf --- /dev/null +++ b/test/prune_old_LA.jl @@ -0,0 +1,69 @@ +methods_to_delete = +[ +:adjoint +:transpose +:inv +:literal_pow +:\ +:/ +:isapprox +:copyto! +:* +:muladd +:copyto! +:isone +:kron! +:kron +:^ +:exp +:cis +:log +:sqrt +:cbrt +:inv +:cos +:sin +:sincos +:tan +:cosh +:sinh +:tanh +:acos +:asin +:atan +:acosh +:asinh +:atanh +:sec +:sech +:csc +:csch +:cot +:coth +:asec +:asech +:acsc +:acot +:acoth +:acsch +] + +let + LA = get(Base.loaded_modules, Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra"), nothing) + if LA !== nothing + @assert hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}}) + for methss in methods_to_delete + meths = getglobal(Base, methss) + for meth in methods(meths) + if meth.module === LA + Base.delete_method(meth) + end + end + end + end + Base.unreference_module(Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra")) +end + +@assert !hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}}) + +pruned_old_LA = true diff --git a/test/qr.jl b/test/qr.jl index b6e9ce3a..2b121ca4 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -2,6 +2,8 @@ module TestQR +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, rmul!, lmul! diff --git a/test/runtests.jl b/test/runtests.jl index a64884dd..0466b047 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,7 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license + +include("prune_old_LA.jl") + using Test, LinearAlgebra for file in readlines(joinpath(@__DIR__, "testgroups")) diff --git a/test/schur.jl b/test/schur.jl index f3d494fb..e0a16b80 100644 --- a/test/schur.jl +++ b/test/schur.jl @@ -2,6 +2,8 @@ module TestSchur +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted diff --git a/test/special.jl b/test/special.jl index ac76279f..3e7a50c5 100644 --- a/test/special.jl +++ b/test/special.jl @@ -2,6 +2,8 @@ module TestSpecial +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: rmul!, BandIndex diff --git a/test/structuredbroadcast.jl b/test/structuredbroadcast.jl index 641c761d..6521ed58 100644 --- a/test/structuredbroadcast.jl +++ b/test/structuredbroadcast.jl @@ -1,6 +1,9 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license module TestStructuredBroadcast + +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") diff --git a/test/svd.jl b/test/svd.jl index 33c02517..fdcfcb30 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -2,6 +2,8 @@ module TestSVD +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted diff --git a/test/symmetric.jl b/test/symmetric.jl index 68f0ff47..e119eba7 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -2,6 +2,8 @@ module TestSymmetric +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") diff --git a/test/symmetriceigen.jl b/test/symmetriceigen.jl index 7a2b88a5..36b69d29 100644 --- a/test/symmetriceigen.jl +++ b/test/symmetriceigen.jl @@ -2,6 +2,8 @@ module TestSymmetricEigen +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations diff --git a/test/triangular.jl b/test/triangular.jl index d63f8b06..a44d3b48 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -2,6 +2,8 @@ module TestTriangular +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: errorbounds, transpose!, BandIndex diff --git a/test/triangular2.jl b/test/triangular2.jl index 4cad1373..f7e761af 100644 --- a/test/triangular2.jl +++ b/test/triangular2.jl @@ -2,6 +2,8 @@ module TestTriangularReal +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasFloat, errorbounds, full!, transpose! diff --git a/test/triangular3.jl b/test/triangular3.jl index 82a2a147..a1990b57 100644 --- a/test/triangular3.jl +++ b/test/triangular3.jl @@ -2,6 +2,8 @@ module TestTriangularComplex +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random using LinearAlgebra: BlasFloat, errorbounds, full!, transpose! diff --git a/test/tridiag.jl b/test/tridiag.jl index 4b592a87..a7f0d7ad 100644 --- a/test/tridiag.jl +++ b/test/tridiag.jl @@ -2,6 +2,8 @@ module TestTridiagonal +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") diff --git a/test/uniformscaling.jl b/test/uniformscaling.jl index 10d427d1..21bd4b80 100644 --- a/test/uniformscaling.jl +++ b/test/uniformscaling.jl @@ -2,6 +2,8 @@ module TestUniformscaling +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") diff --git a/test/unitful.jl b/test/unitful.jl index cd4ade52..84e875e7 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -1,5 +1,7 @@ module TestUnitfulLinAlg +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra, Random Random.seed!(1234321) From 14c6ad37243df443fde8180abd87950b842085c1 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 1 May 2025 17:18:37 +0200 Subject: [PATCH 18/28] Fix multiplication with empty `HessenbergQ` (#1326) (cherry picked from commit b14390e2aec1d4fe2ede9ba2992d104d2d3759e2) --- src/lapack.jl | 5 +++-- test/hessenberg.jl | 7 +++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/lapack.jl b/src/lapack.jl index a0c1c785..20af7ca8 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -6130,10 +6130,11 @@ for (ormhr, elty) in chkside(side) chktrans(trans) n = checksquare(A) + ntau = length(tau) mC, nC = size(C, 1), size(C, 2) - if n - length(tau) != 1 - throw(DimensionMismatch(lazy"tau has length $(length(tau)), needs $(n - 1)")) + if max(n - 1, 0) != ntau + throw(DimensionMismatch(lazy"tau has length $ntau, needs $(n - 1)")) end if (side == 'L' && mC != n) || (side == 'R' && nC != n) throw(DimensionMismatch("A and C matrices are not conformable")) diff --git a/test/hessenberg.jl b/test/hessenberg.jl index 91f6d155..8485c014 100644 --- a/test/hessenberg.jl +++ b/test/hessenberg.jl @@ -286,4 +286,11 @@ end end end +@testset "multiplication with empty HessenbergQ" begin + @test ones(2, 0)*hessenberg(zeros(0,0)).Q == zeros(2,0) + @test_throws DimensionMismatch ones(2, 1)*hessenberg(zeros(0,0)).Q + @test hessenberg(zeros(0,0)).Q * ones(0, 2) == zeros(0,2) + @test_throws DimensionMismatch hessenberg(zeros(0,0)).Q * ones(1, 2) +end + end # module TestHessenberg From 380f0668ccd637466c0bbd123b497061d4546588 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 11 May 2025 21:37:19 +0530 Subject: [PATCH 19/28] Test: prune old LA based on ENV variable (#1335) This PR prunes the old `LinearAlgebra` module based on an environment variable that is set in `.ci/run_tests.jl`. Hopefully, this would mean that tests that are being run without this environment variable would run as expected without precompilation issues. The goal of this PR is to ensure that when julia is being built with the master branch of `LinearAlgebra`, the tests would work without the need for pruning (i.e., we want to fix the test failures in https://github.com/JuliaLang/julia/pull/58242). This would also mean that to run the tests locally, one would need to set the `JULIA_PRUNE_OLD_LA` environment variable to `true`. (cherry picked from commit 87d4c8b08e456eaf88a088945576e30c7660350d) --- .ci/run_tests.jl | 5 +- README.md | 31 +++++++++++- docs/make.jl | 4 +- test/bitarray.jl | 6 +++ test/prune_old_LA.jl | 118 +++++++++++++++++++++++-------------------- 5 files changed, 104 insertions(+), 60 deletions(-) diff --git a/.ci/run_tests.jl b/.ci/run_tests.jl index 8a54ef62..879b73b1 100644 --- a/.ci/run_tests.jl +++ b/.ci/run_tests.jl @@ -6,6 +6,5 @@ end proj = abspath(joinpath(@__DIR__, "..")) cmd = """Base.runtests(["LinearAlgebra"]; propagate_project=true, ncores=$ncores)""" -withenv("JULIA_NUM_THREADS" => 1) do - run(`$(Base.julia_cmd()) --project=$proj --compiled-modules=existing -e $cmd`) -end +run(addenv(`$(Base.julia_cmd()) --project=$proj --compiled-modules=existing -e $cmd`, + "JULIA_NUM_THREADS" => 1, "JULIA_PRUNE_OLD_LA" => true)) diff --git a/README.md b/README.md index a39a7248..c0951dcc 100644 --- a/README.md +++ b/README.md @@ -37,10 +37,37 @@ This package performs some type piracy and is also included in the sysimage, whi To use a development version of this package, you can choose one of the following methods: 1. **Change the UUID in the project file and load the package:** - This approach will produce warnings and may lead to method ambiguities between the development version and the one in the sysimage, but it can be used for basic experimentation. + + Change the UUID line in `Project.toml` as + ```diff + - uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + + uuid = "27e2e46d-f89d-539d-b4ee-838fcccc9c8e" + ``` + + Start `julia` as + ```console + JULIA_PRUNE_OLD_LA=true julia +nightly --compiled-modules=existing --project + ``` + where it is assumed that one is already within the `LinearAlgebra` directory (otherwise, adjust + the project path accordingly). The `julia +nightly` command above assumes that `juliaup` is being used + to launch `julia`, but one may substitute this with the path to the julia executable. + + Within the `julia` session, load the `LinearAlgebra` module after pruning the one in the sysimage. This may be done as + ```julia + include("test/prune_old_LA.jl") && using LinearAlgebra + ``` + + Note that loading the test files in the REPL will automatically carry out the pruning to ensure that the development version of the package is used in the tests. + + If you are contributing to the repo using this method, it may be convenient to ignore the local changes to `Project.toml` by running + ```console + git update-index --skip-worktree Project.toml + ``` 2. **Build Julia with the custom `LinearAlgebra` commit:** - Modify the commit in `stdlib/LinearAlgebra.version` and build Julia. + + Modify the commit in [`stdlib/LinearAlgebra.version`](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra.version) and build Julia. This requires one to push the development branch + to `GitHub` or an equivalent platform. 3. **Build a custom sysimage with the new `LinearAlgebra`:** - Install `PackageCompiler`. diff --git a/docs/make.jl b/docs/make.jl index a1d3a2d0..0f8b79d9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,4 +1,6 @@ -include("../test/prune_old_LA.jl") +withenv("JULIA_PRUNE_OLD_LA" => "true") do + include("../test/prune_old_LA.jl") +end using LinearAlgebra using Documenter: DocMeta, makedocs, deploydocs, HTML diff --git a/test/bitarray.jl b/test/bitarray.jl index a185873e..10b9d4b8 100644 --- a/test/bitarray.jl +++ b/test/bitarray.jl @@ -1,3 +1,7 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +module TestBitArray + isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") using LinearAlgebra, Test, Random @@ -95,3 +99,5 @@ b2 = bitrand(v1) b1 = bitrand(n1, n1) @check_bit_operation diag(b1) + +end # module diff --git a/test/prune_old_LA.jl b/test/prune_old_LA.jl index f628cddf..7e64cfa3 100644 --- a/test/prune_old_LA.jl +++ b/test/prune_old_LA.jl @@ -1,56 +1,58 @@ -methods_to_delete = -[ -:adjoint -:transpose -:inv -:literal_pow -:\ -:/ -:isapprox -:copyto! -:* -:muladd -:copyto! -:isone -:kron! -:kron -:^ -:exp -:cis -:log -:sqrt -:cbrt -:inv -:cos -:sin -:sincos -:tan -:cosh -:sinh -:tanh -:acos -:asin -:atan -:acosh -:asinh -:atanh -:sec -:sech -:csc -:csch -:cot -:coth -:asec -:asech -:acsc -:acot -:acoth -:acsch -] - let - LA = get(Base.loaded_modules, Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra"), nothing) - if LA !== nothing + methods_to_delete = + [ + :adjoint + :transpose + :inv + :literal_pow + :\ + :/ + :isapprox + :copyto! + :* + :muladd + :copyto! + :isone + :kron! + :kron + :^ + :exp + :cis + :log + :sqrt + :cbrt + :inv + :cos + :sin + :sincos + :tan + :cosh + :sinh + :tanh + :acos + :asin + :atan + :acosh + :asinh + :atanh + :sec + :sech + :csc + :csch + :cot + :coth + :asec + :asech + :acsc + :acot + :acoth + :acsch + ] + + prune_old_LA = parse(Bool, get(ENV, "JULIA_PRUNE_OLD_LA", "false")) + LinalgSysImg = Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra") + LA = get(Base.loaded_modules, LinalgSysImg, nothing) + if LA !== nothing && prune_old_LA @assert hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}}) for methss in methods_to_delete meths = getglobal(Base, methss) @@ -61,9 +63,17 @@ let end end end - Base.unreference_module(Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra")) end -@assert !hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}}) +# check in a separate block to ensure that the latest world age is used +let + prune_old_LA = parse(Bool, get(ENV, "JULIA_PRUNE_OLD_LA", "false")) + LinalgSysImg = Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra") + LA = get(Base.loaded_modules, LinalgSysImg, nothing) + if LA !== nothing && prune_old_LA + @assert !hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}}) + end + prune_old_LA && Base.unreference_module(LinalgSysImg) +end pruned_old_LA = true From fda7486c0ac30ea49bda9271bd67c26f07b0c8eb Mon Sep 17 00:00:00 2001 From: Michael Abbott <32575566+mcabbott@users.noreply.github.com> Date: Tue, 13 May 2025 00:27:32 -0400 Subject: [PATCH 20/28] Add compat notice for `diagview` (#1355) This function didn't exist in Julia 1.11 There is a note in NEWS, here: https://github.com/JuliaLang/julia/blob/master/HISTORY.md#linearalgebra Added in https://github.com/JuliaLang/julia/pull/56175 (cherry picked from commit 2cb1e9862922630628054049381a4eb692fb0eaa) --- src/dense.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/dense.jl b/src/dense.jl index 234391e5..84941d3e 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -297,6 +297,9 @@ Return a view into the `k`th diagonal of the matrix `M`. See also [`diag`](@ref), [`diagind`](@ref). +!!! compat "Julia 1.12" + This function requires Julia 1.12 or later. + # Examples ```jldoctest julia> A = [1 2 3; 4 5 6; 7 8 9] From b910d5ecb9ca4508e82a2b21d30c627dca52438e Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 May 2025 12:26:03 +0530 Subject: [PATCH 21/28] Prune `LinearAlgebra` module in ambiguity test (#1349) The ambiguity check was being run in a separate process, but it was not pruning the `LinearAlgebra` module. This was therefore not testing for ambiguities in the latest code, but the one in the sysimg. The PR would be able to replicate the test failure seen in https://buildkite.com/julialang/julia-master/builds/47424#0196c03a-2230-4d64-b404-e04dbfc76c3f (cherry picked from commit 0b46d5c8c9b53ff8b0fcbf7e23a8f254ff34fc15) --- test/ambiguous_exec.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/ambiguous_exec.jl b/test/ambiguous_exec.jl index 7b89c0a4..09693791 100644 --- a/test/ambiguous_exec.jl +++ b/test/ambiguous_exec.jl @@ -1,5 +1,9 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license +module TestAmbiguity + +isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") + using Test, LinearAlgebra let ambig = detect_ambiguities(LinearAlgebra; recursive=true) @test isempty(ambig) @@ -19,3 +23,5 @@ let ambig = detect_ambiguities(LinearAlgebra; recursive=true) @test isempty(expect) @test good end + +end # module From 465b5b915de2cad7555aee9a472f997520695199 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 May 2025 12:54:58 +0530 Subject: [PATCH 22/28] Ensure positive-definite matrix in lapack posv test (#1238) --- test/lapack.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/lapack.jl b/test/lapack.jl index 6e7bab0b..fdf5b6bd 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -679,7 +679,7 @@ end @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) local n = 10 a = randn(elty, n, n) - A = a'*a + A = a'*a + I B = rand(elty, n, n) D = copy(A) C = copy(B) From 69fa500562a856fef5062adb4722844b591318cc Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 May 2025 13:00:25 +0530 Subject: [PATCH 23/28] Add diagm example (#1298) --- src/dense.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/dense.jl b/src/dense.jl index 84941d3e..636044d3 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -402,6 +402,7 @@ Construct a matrix with elements of the vector as diagonal elements. By default, the matrix is square and its size is given by `length(v)`, but a non-square size `m`×`n` can be specified by passing `m,n` as the first arguments. +The diagonal will be zero-padded if necessary. # Examples ```jldoctest @@ -410,6 +411,13 @@ julia> diagm([1,2,3]) 1 0 0 0 2 0 0 0 3 + +julia> diagm(4, 5, [1,2,3]) +4×5 Matrix{Int64}: + 1 0 0 0 0 + 0 2 0 0 0 + 0 0 3 0 0 + 0 0 0 0 0 ``` """ diagm(v::AbstractVector) = diagm(0 => v) From f6d16c220433665991b78e2623a440a9a9906282 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 28 Apr 2025 23:45:36 +0530 Subject: [PATCH 24/28] Copy matrices in `triu`/`tril` if no zero exists for the `eltype` (#1320) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This fixes a regression in the `triu`/`tril` implementation on v1.12 and nightly. After this, the following works again: ```julia julia> M = fill(ones(2,2), 3, 3) 3×3 Matrix{Matrix{Float64}}: [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] julia> using Test: GenericArray julia> triu(GenericArray(M),1) 3×3 GenericArray{Matrix{Float64}, 2}: [0.0 0.0; 0.0 0.0] [1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0] [0.0 0.0; 0.0 0.0] [0.0 0.0; 0.0 0.0] [1.0 1.0; 1.0 1.0] [0.0 0.0; 0.0 0.0] [0.0 0.0; 0.0 0.0] [0.0 0.0; 0.0 0.0] ``` --- src/generic.jl | 22 ++++++++++++++++++++-- test/dense.jl | 10 ++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index fc0b4664..77740d9f 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -446,7 +446,8 @@ julia> triu(a,-3) 1.0 1.0 1.0 1.0 ``` """ -function triu(M::AbstractMatrix, k::Integer = 0) +triu(M::AbstractMatrix, k::Integer = 0) = _triu(M, Val(haszero(eltype(M))), k) +function _triu(M::AbstractMatrix, ::Val{true}, k::Integer) d = similar(M) A = triu!(d,k) if iszero(k) @@ -459,6 +460,14 @@ function triu(M::AbstractMatrix, k::Integer = 0) end return A end +function _triu(M::AbstractMatrix, ::Val{false}, k::Integer) + d = similar(M) + # since the zero would need to be evaluated from the elements, + # we copy the array to avoid undefined references in triu! + copy!(d, M) + A = triu!(d,k) + return A +end """ tril(M, k::Integer = 0) @@ -489,7 +498,8 @@ julia> tril(a,-3) 1.0 0.0 0.0 0.0 ``` """ -function tril(M::AbstractMatrix,k::Integer=0) +tril(M::AbstractMatrix,k::Integer=0) = _tril(M, Val(haszero(eltype(M))), k) +function _tril(M::AbstractMatrix, ::Val{true}, k::Integer) d = similar(M) A = tril!(d,k) if iszero(k) @@ -502,6 +512,14 @@ function tril(M::AbstractMatrix,k::Integer=0) end return A end +function _tril(M::AbstractMatrix, ::Val{false}, k::Integer) + d = similar(M) + # since the zero would need to be evaluated from the elements, + # we copy the array to avoid undefined references in tril! + copy!(d, M) + A = tril!(d,k) + return A +end """ triu!(M) diff --git a/test/dense.jl b/test/dense.jl index 47efb78c..743c0ed1 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1374,4 +1374,14 @@ end @test !any(LinearAlgebra.getstructure(A)) end +@testset "triu/tril for block matrices" begin + O = ones(2,2) + Z = zero(O) + M = fill(O, 3, 3) + res = fill(Z, size(M)) + res[1,2] = res[1,3] = res[2,3] = O + @test triu(GenericArray(M),1) == res + @test tril(GenericArray(M),-1) == res' +end + end # module TestDense From f1378583686c9c54d1a32b81c195e50e72dea16a Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 6 May 2025 08:51:10 +0530 Subject: [PATCH 25/28] Make `fillstored!` public (#1333) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `fillstored!` offers a way to fill only the stored indices of a structured matrix without having to populate the parent. E.g.: ```julia julia> U = UpperTriangular(zeros(4,4)) 4×4 UpperTriangular{Float64, Matrix{Float64}}: 0.0 0.0 0.0 0.0 ⋅ 0.0 0.0 0.0 ⋅ ⋅ 0.0 0.0 ⋅ ⋅ ⋅ 0.0 julia> LinearAlgebra.fillstored!(U, 2) 4×4 UpperTriangular{Float64, Matrix{Float64}}: 2.0 2.0 2.0 2.0 ⋅ 2.0 2.0 2.0 ⋅ ⋅ 2.0 2.0 ⋅ ⋅ ⋅ 2.0 ``` This seems like a useful function that should be public. This came up on discourse: https://discourse.julialang.org/t/how-to-set-all-elements-in-a-lower-triangular-matrix/128603 --- src/LinearAlgebra.jl | 3 ++- src/special.jl | 25 ++++++++++++++++++++++++- src/symmetric.jl | 10 ++++------ test/special.jl | 24 ++++++++++++++++++++++++ test/symmetric.jl | 9 +++++++++ 5 files changed, 63 insertions(+), 8 deletions(-) diff --git a/src/LinearAlgebra.jl b/src/LinearAlgebra.jl index d2f82705..6f6e21ba 100644 --- a/src/LinearAlgebra.jl +++ b/src/LinearAlgebra.jl @@ -179,7 +179,8 @@ public AbstractTriangular, symmetric, symmetric_type, zeroslike, - matprod_dest + matprod_dest, + fillstored! const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32} const BlasReal = Union{Float64,Float32} diff --git a/src/special.jl b/src/special.jl index 58863858..1327af67 100644 --- a/src/special.jl +++ b/src/special.jl @@ -320,10 +320,33 @@ _diag_or_value(A::Diagonal) = A.diag _diag_or_value(A::UniformScaling) = A.λ # fill[stored]! methods +""" + fillstored!(A::AbstractMatrix, x) + +Fill only the stored indices of a structured matrix `A` with the value `x`. + +# Example +```jldoctest +julia> A = Tridiagonal(zeros(2), zeros(3), zeros(2)) +3×3 Tridiagonal{Float64, Vector{Float64}}: + 0.0 0.0 ⋅ + 0.0 0.0 0.0 + ⋅ 0.0 0.0 + +julia> LinearAlgebra.fillstored!(A, 2) +3×3 Tridiagonal{Float64, Vector{Float64}}: + 2.0 2.0 ⋅ + 2.0 2.0 2.0 + ⋅ 2.0 2.0 +``` +""" fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A) fillstored!(A::Bidiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A) fillstored!(A::Tridiagonal, x) = (fill!(A.dl, x); fill!(A.d, x); fill!(A.du, x); A) -fillstored!(A::SymTridiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A) +function fillstored!(A::SymTridiagonal, x) + issymmetric(x) || throw(ArgumentError("cannot set a diagonal entry of a SymTridiagonal to an asymmetric value")) + (fill!(A.dv, x); fill!(A.ev, x); A) +end _small_enough(A::Union{Diagonal, Bidiagonal}) = size(A, 1) <= 1 _small_enough(A::Tridiagonal) = size(A, 1) <= 2 diff --git a/src/symmetric.jl b/src/symmetric.jl index 11b34a99..1edb0570 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -401,13 +401,11 @@ fill!(A::HermOrSym, x) = fillstored!(A, x) function fillstored!(A::HermOrSym{T}, x) where T xT = convert(T, x) if isa(A, Hermitian) - isreal(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a nonreal value")) - end - if A.uplo == 'U' - fillband!(A.data, xT, 0, size(A,2)-1) - else # A.uplo == 'L' - fillband!(A.data, xT, 1-size(A,1), 0) + ishermitian(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a non-hermitian value")) + elseif isa(A, Symmetric) + issymmetric(xT) || throw(ArgumentError("cannot fill Symmetric matrix with an asymmetric value")) end + applytri(A -> fillstored!(A, xT), A) return A end diff --git a/test/special.jl b/test/special.jl index 3e7a50c5..0cb85d56 100644 --- a/test/special.jl +++ b/test/special.jl @@ -835,4 +835,28 @@ end end end +@testset "fillstored!" begin + dv, ev = zeros(4), zeros(3) + D = Diagonal(dv) + LinearAlgebra.fillstored!(D, 2) + @test D == diagm(fill(2, length(dv))) + + dv .= 0 + B = Bidiagonal(dv, ev, :U) + LinearAlgebra.fillstored!(B, 2) + @test B == diagm(0=>fill(2, length(dv)), 1=>fill(2, length(ev))) + + dv .= 0 + ev .= 0 + T = Tridiagonal(ev, dv, ev) + LinearAlgebra.fillstored!(T, 2) + @test T == diagm(-1=>fill(2, length(ev)), 0=>fill(2, length(dv)), 1=>fill(2, length(ev))) + + dv .= 0 + ev .= 0 + ST = SymTridiagonal(dv, ev) + LinearAlgebra.fillstored!(ST, 2) + @test ST == diagm(-1=>fill(2, length(ev)), 0=>fill(2, length(dv)), 1=>fill(2, length(ev))) +end + end # module TestSpecial diff --git a/test/symmetric.jl b/test/symmetric.jl index e119eba7..17500537 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -1165,4 +1165,13 @@ end end end +@testset "fillstored!" begin + A = zeros(4,4) + for T in (Symmetric, Hermitian) + A .= 0 + LinearAlgebra.fillstored!(T(A), 2) + @test all(==(2), T(A)) + end +end + end # module TestSymmetric From aefab75e8d4de0326b98ee29c9acacce1d2f462b Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 May 2025 13:16:09 +0530 Subject: [PATCH 26/28] Document SingularException throw for inv(::AbstractMatrix) (#1331) --- src/generic.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/generic.jl b/src/generic.jl index 77740d9f..7c69f20c 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1134,6 +1134,8 @@ Matrix inverse. Computes matrix `N` such that Computed by solving the left-division `N = M \\ I`. +A [`SingularException`](@ref) is thrown if `M` fails numerical inversion. + # Examples ```jldoctest julia> M = [2 5; 1 3] From e6e7449c7aa5cc5739decdb453b930871bbb4efd Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 May 2025 18:55:36 +0530 Subject: [PATCH 27/28] Fix copy for partly initialized unit triangular (#1350) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Currently, we forward copy for unit triangular matrices to the parents. This also copies the diagonal, which may not be initialized in the parent. This PR fixes this. After this, the following works again ```julia julia> M = Matrix{BigFloat}(undef,2,2); julia> M[2,1] = 3; julia> M' 2×2 adjoint(::Matrix{BigFloat}) with eltype BigFloat: #undef 3.0 #undef #undef julia> U = UnitUpperTriangular(M') 2×2 UnitUpperTriangular{BigFloat, Adjoint{BigFloat, Matrix{BigFloat}}}: 1.0 3.0 ⋅ 1.0 julia> copyto!(similar(M), U) 2×2 Matrix{BigFloat}: 1.0 3.0 0.0 1.0 ``` --- src/triangular.jl | 7 ++++--- test/triangular.jl | 11 +++++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/triangular.jl b/src/triangular.jl index c4da1adb..6ce33dc7 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -599,14 +599,15 @@ end copytrito!(dest, U, U isa UpperOrUnitUpperTriangular ? 'L' : 'U') return dest end -@propagate_inbounds function _copyto!(dest::StridedMatrix, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}) +@propagate_inbounds function _copyto!(dest::StridedMatrix, + U::Union{UnitUpperOrUnitLowerTriangular, UpperOrLowerTriangular{<:Any, <:StridedMatrix}}) U2 = Base.unalias(dest, U) copyto_unaliased!(dest, U2) return dest end # for strided matrices, we explicitly loop over the arrays to improve cache locality # This fuses the copytrito! for the two halves -@inline function copyto_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular{<:Any, <:StridedMatrix}) +@inline function copyto_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular) @boundscheck checkbounds(dest, axes(U)...) isunit = U isa UnitUpperTriangular for col in axes(dest,2) @@ -619,7 +620,7 @@ end end return dest end -@inline function copyto_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular{<:Any, <:StridedMatrix}) +@inline function copyto_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular) @boundscheck checkbounds(dest, axes(L)...) isunit = L isa UnitLowerTriangular for col in axes(dest,2) diff --git a/test/triangular.jl b/test/triangular.jl index a44d3b48..1a713747 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -6,6 +6,7 @@ isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl") using Test, LinearAlgebra, Random using LinearAlgebra: errorbounds, transpose!, BandIndex +using Test: GenericArray const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") @@ -650,6 +651,16 @@ end @test_throws "cannot set index in the upper triangular part" copyto!(A, B) end end + + @testset "partly initialized unit triangular" begin + for T in (UnitUpperTriangular, UnitLowerTriangular) + isupper = T == UnitUpperTriangular + M = Matrix{BigFloat}(undef, 2, 2) + M[1+!isupper,1+isupper] = 3 + U = T(GenericArray(M)) + @test copyto!(similar(M), U) == U + end + end end @testset "getindex with Integers" begin From edca55925a4bdd3849d9c71f4ca29c9513af632f Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 16 May 2025 22:50:58 +0530 Subject: [PATCH 28/28] log for dense diagonal matrix with negative elements (#1352) Fixes a regression introduced in v1.12, where we have a special branch for diagonal matrices to improve performance. The issue is that log for `Diagonal` matrices acts elementwise, which, for real matrices, requires the diagonal elements to be greater than zero. This PR adds an extra branch to check this. --- src/dense.jl | 8 ++++++-- test/dense.jl | 11 +++++++++++ 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/src/dense.jl b/src/dense.jl index 636044d3..234b1022 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -901,8 +901,12 @@ julia> log(A) """ function log(A::AbstractMatrix) # If possible, use diagonalization - if isdiag(A) - return applydiagonal(log, A) + if isdiag(A) && eltype(A) <: Union{Real,Complex} + if eltype(A) <: Real && all(>=(0), diagview(A)) + return applydiagonal(log, A) + else + return applydiagonal(log∘complex, A) + end elseif ishermitian(A) logHermA = log(Hermitian(A)) return ishermitian(logHermA) ? copytri!(parent(logHermA), 'U', true) : parent(logHermA) diff --git a/test/dense.jl b/test/dense.jl index 743c0ed1..d0cfa007 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1384,4 +1384,15 @@ end @test tril(GenericArray(M),-1) == res' end +@testset "log for diagonal" begin + D = diagm([-2.0, 2.0]) + @test log(D) ≈ log(UpperTriangular(D)) + D = diagm([-2.0, 0.0]) + @test log(D) ≈ log(UpperTriangular(D)) + D = diagm([2.0, 2.0]) + @test log(D) ≈ log(UpperTriangular(D)) + D = diagm([2.0, 2.0*im]) + @test log(D) ≈ log(UpperTriangular(D)) +end + end # module TestDense