From 9941132efd68eb5d5df5d005a4b4eb73b02135c7 Mon Sep 17 00:00:00 2001 From: Luis-Varona Date: Wed, 4 Jun 2025 10:35:05 -0300 Subject: [PATCH 01/12] Add docstring to the rank(::QRPivoted) method With the release of Julia 1.12, LinearAlgebra.jl is adding dispatches for rank to SVD and QRPivoted objects, allowing the re-use of existing matrix factorizations where they have already been computed. The SVD method already has an appropriate docstring; this PR adds a near-identical one to the QRPivoted method as well, with some additional context for how it differs from the default SVD-based computation. --- src/qr.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/qr.jl b/src/qr.jl index 9e0f49d1..5c59ab60 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -535,6 +535,28 @@ function ldiv!(A::QRCompactWY{T}, B::AbstractMatrix{T}) where {T} return B end +""" + rank(A::QRPivoted{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T} + +Compute the numerical rank of the QR factorization `A` by counting how many diagonal entries of +`A.factors` are greater than `max(atol, rtol*Δ₁)` where `Δ₁` is the largest calculated such entry. +This is similar to the [`rank(::AbstractMatrix)`](@ref) method insofar as it counts the number of +(numerically) nonzero coefficients from a matrix factorization, although the default method uses an +SVD instead of a QR factorization. Like [`rank(::SVD)`](@ref), this method also re-uses an existing +matrix factorization. + +Using a QR factorization to compute rank should typically produce the same result as using SVD, +although it may be more prone to overestimating the rank in pathological cases where the matrix is +ill-conditioned. It is also worth noting that it is generally faster to compute a QR factorization +than it is to compute an SVD, so this method may be preferred when performance is a concern. + +`atol` and `rtol` are the absolute and relative tolerances, respectively. +The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `A` +and `ϵ` is the [`eps`](@ref) of the element type of `A`. + +!!! compat "Julia 1.12" + The `rank(::QRPivoted)` method requires at least Julia 1.12. +""" function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(eltype(A)))) * iszero(atol)) m = min(size(A)...) m == 0 && return 0 From 213b31d165ec31d85633690faba4bfc52f26e75a Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Wed, 13 Aug 2025 09:11:51 -0400 Subject: [PATCH 02/12] improve type stability of structural zero error for triangular (#1413) Co-authored-by: Viral B. Shah Co-authored-by: Jishnu Bhattacharya --- src/triangular.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/triangular.jl b/src/triangular.jl index 6ce33dc7..70e7b58d 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -248,16 +248,15 @@ Base.@constprop :aggressive @propagate_inbounds function getindex(A::Union{Lower _shouldforwardindex(A, b) ? A.data[b] : diagzero(A.data, b) end -_zero_triangular_half_str(::Type{<:UpperOrUnitUpperTriangular}) = "lower" -_zero_triangular_half_str(::Type{<:LowerOrUnitLowerTriangular}) = "upper" +_zero_triangular_half_str(T::Type) = T <: UpperOrUnitUpperTriangular ? "lower" : "upper" -@noinline function throw_nonzeroerror(T, @nospecialize(x), i, j) +@noinline function throw_nonzeroerror(T::DataType, @nospecialize(x), i, j) Ts = _zero_triangular_half_str(T) Tn = nameof(T) throw(ArgumentError( lazy"cannot set index in the $Ts triangular part ($i, $j) of an $Tn matrix to a nonzero value ($x)")) end -@noinline function throw_nononeerror(T, @nospecialize(x), i, j) +@noinline function throw_nononeerror(T::DataType, @nospecialize(x), i, j) Tn = nameof(T) throw(ArgumentError( lazy"cannot set index on the diagonal ($i, $j) of an $Tn matrix to a non-unit value ($x)")) @@ -303,7 +302,7 @@ end return A end -@noinline function throw_setindex_structuralzero_error(T, @nospecialize(x)) +@noinline function throw_setindex_structuralzero_error(T::DataType, @nospecialize(x)) Ts = _zero_triangular_half_str(T) Tn = nameof(T) throw(ArgumentError( From f26ced15ab6e65a722edbe7c3724c38a624901c1 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 16 Jun 2025 12:27:15 +0530 Subject: [PATCH 03/12] Forward `diagview` for `adjoint`/`transpose` to parent (#1373) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit After this, ```julia julia> B = reshape([1:9;], 3, 3) 3×3 Matrix{Int64}: 1 4 7 2 5 8 3 6 9 julia> diagview(B', 1) 2-element view(::Vector{Int64}, 2:4:6) with eltype Int64: 2 6 ``` The view directly indexes into the parent, and the `adjoint` layer is removed. This also makes the view linearly indexed. --- src/adjtrans.jl | 20 ++++++++++++++++++++ test/adjtrans.jl | 11 +++++++++++ 2 files changed, 31 insertions(+) diff --git a/src/adjtrans.jl b/src/adjtrans.jl index d81aa3ae..c041d50d 100644 --- a/src/adjtrans.jl +++ b/src/adjtrans.jl @@ -536,3 +536,23 @@ conj(A::Adjoint) = transpose(A.parent) function Base.replace_in_print_matrix(A::AdjOrTrans,i::Integer,j::Integer,s::AbstractString) Base.replace_in_print_matrix(parent(A), j, i, s) end + +# Special adjoint/transpose methods for Adjoint/Transpose that have been reshaped as a vector +# these are used in transposing banded matrices by forwarding the operation to the bands +""" + _vectranspose(A::AbstractVector)::AbstractVector + +Compute `vec(transpose(A))`, but avoid an allocating reshape if possible +""" +_vectranspose(A::AbstractVector) = vec(transpose(A)) +_vectranspose(A::Base.ReshapedArray{<:Any,1,<:TransposeAbsVec}) = transpose(parent(A)) +""" + _vecadjoint(A::AbstractVector)::AbstractVector + +Compute `vec(adjoint(A))`, but avoid an allocating reshape if possible +""" +_vecadjoint(A::AbstractVector) = vec(adjoint(A)) +_vecadjoint(A::Base.ReshapedArray{<:Any,1,<:AdjointAbsVec}) = adjoint(parent(A)) + +diagview(A::Transpose, k::Integer = 0) = _vectranspose(diagview(parent(A), -k)) +diagview(A::Adjoint, k::Integer = 0) = _vecadjoint(diagview(parent(A), -k)) diff --git a/test/adjtrans.jl b/test/adjtrans.jl index bdeaa697..4de72fd6 100644 --- a/test/adjtrans.jl +++ b/test/adjtrans.jl @@ -749,4 +749,15 @@ end end end +@testset "diagview" begin + for A in (rand(4, 4), rand(ComplexF64,4,4), + fill([1 2; 3 4], 4, 4)) + for k in -3:3 + @test diagview(A', k) == diag(A', k) + @test diagview(transpose(A), k) == diag(transpose(A), k) + end + @test IndexStyle(diagview(A')) == IndexLinear() + end +end + end # module TestAdjointTranspose From 29a195cf2d5a0e56d779c563fe1ad803ebba92f1 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 18 Jun 2025 13:43:23 +0530 Subject: [PATCH 04/12] Specialize `triu!` and `tril!` for adj/trans (#1385) We may forward the operation to the parent, which would usually ensure that the loops are cache-friendly. --- src/adjtrans.jl | 4 ++++ test/adjtrans.jl | 11 +++++++++++ 2 files changed, 15 insertions(+) diff --git a/src/adjtrans.jl b/src/adjtrans.jl index c041d50d..ee610be7 100644 --- a/src/adjtrans.jl +++ b/src/adjtrans.jl @@ -556,3 +556,7 @@ _vecadjoint(A::Base.ReshapedArray{<:Any,1,<:AdjointAbsVec}) = adjoint(parent(A)) diagview(A::Transpose, k::Integer = 0) = _vectranspose(diagview(parent(A), -k)) diagview(A::Adjoint, k::Integer = 0) = _vecadjoint(diagview(parent(A), -k)) + +# triu and tril +triu!(A::AdjOrTransAbsMat, k::Integer = 0) = wrapperop(A)(tril!(parent(A), -k)) +tril!(A::AdjOrTransAbsMat, k::Integer = 0) = wrapperop(A)(triu!(parent(A), -k)) diff --git a/test/adjtrans.jl b/test/adjtrans.jl index 4de72fd6..e5654004 100644 --- a/test/adjtrans.jl +++ b/test/adjtrans.jl @@ -760,4 +760,15 @@ end end end +@testset "triu!/tril!" begin + @testset for sz in ((4,4), (3,4), (4,3)) + A = rand(sz...) + B = similar(A) + @testset for f in (adjoint, transpose), k in -3:3 + @test triu!(f(copy!(B, A)), k) == triu(f(A), k) + @test tril!(f(copy!(B, A)), k) == tril!(f(A), k) + end + end +end + end # module TestAdjointTranspose From fb13e94d527351510b0660ff33041951f4501d61 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 4 Jul 2025 01:31:33 -0400 Subject: [PATCH 05/12] restrict inv(::Factorization) to square matrices (#1397) --- src/factorization.jl | 2 +- src/svd.jl | 1 + test/svd.jl | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 4cefc661..b035f642 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -110,7 +110,7 @@ Factorization{T}(A::AdjointFactorization) where {T} = adjoint(Factorization{T}(parent(A))) Factorization{T}(A::TransposeFactorization) where {T} = transpose(Factorization{T}(parent(A))) -inv(F::Factorization{T}) where {T} = (n = size(F, 1); ldiv!(F, Matrix{T}(I, n, n))) +inv(F::Factorization{T}) where {T} = (n = checksquare(F); ldiv!(F, Matrix{T}(I, n, n))) Base.hash(F::Factorization, h::UInt) = mapreduce(f -> hash(getfield(F, f)), hash, 1:nfields(F); init=h) Base.:(==)( F::T, G::T) where {T<:Factorization} = all(f -> getfield(F, f) == getfield(G, f), 1:nfields(F)) diff --git a/src/svd.jl b/src/svd.jl index 1fa5a719..ce86284b 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -272,6 +272,7 @@ function ldiv!(A::SVD{T}, B::AbstractVecOrMat) where T end function inv(F::SVD{T}) where T + checksquare(F) @inbounds for i in eachindex(F.S) iszero(F.S[i]) && throw(SingularException(i)) end diff --git a/test/svd.jl b/test/svd.jl index fdcfcb30..d73c839a 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -53,9 +53,9 @@ using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted @test sf2.U*Diagonal(sf2.S)*sf2.Vt' ≊ m2 @test ldiv!([0., 0.], svd(Matrix(I, 2, 2)), [1., 1.]) ≊ [1., 1.] + @test_throws DimensionMismatch inv(svd(Matrix(I, 3, 2))) @test inv(svd(Matrix(I, 2, 2))) ≈ I @test inv(svd([1 2; 3 4])) ≈ [-2.0 1.0; 1.5 -0.5] - @test inv(svd([1 0 1; 0 1 0])) ≈ [0.5 0.0; 0.0 1.0; 0.5 0.0] @test_throws SingularException inv(svd([0 0; 0 0])) @test inv(svd([1+2im 3+4im; 5+6im 7+8im])) ≈ [-0.5 + 0.4375im 0.25 - 0.1875im; 0.375 - 0.3125im -0.125 + 0.0625im] end From a23c5304bf4645111aa001bd1c2a724c7ece4577 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Sat, 5 Jul 2025 13:22:03 -0400 Subject: [PATCH 06/12] fix: add BLAS.geru! manual doc (#1402) closes #1401 --- docs/src/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 3e18a457..362f4711 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -745,7 +745,7 @@ LinearAlgebra.BLAS.trsv ```@docs LinearAlgebra.BLAS.ger! -# xGERU +LinearAlgebra.BLAS.geru! # xGERC LinearAlgebra.BLAS.her! # xHPR From 4fe709da8d19ad55a9deba248512ae8c801f242c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 18 Aug 2025 13:34:37 +0200 Subject: [PATCH 07/12] Respect strong zero in generic matvec-mul (#1417) --- src/matmul.jl | 12 +++++++----- test/matmul.jl | 20 ++++++++++++++++---- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/src/matmul.jl b/src/matmul.jl index 10a4c81c..2548a239 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -968,11 +968,13 @@ function __generic_matvecmul!(::typeof(identity), C::AbstractVector, A::Abstract C[i] = zero(A[i]*B[1] + A[i]*B[1]) end end - for k = eachindex(B) - aoffs = (k-1)*Astride - b = @stable_muladdmul MulAddMul(alpha,false)(B[k]) - for i = eachindex(C) - C[i] += A[aoffs + i] * b + if !iszero(alpha) + for k = eachindex(B) + aoffs = (k-1)*Astride + b = @stable_muladdmul MulAddMul(alpha,false)(B[k]) + for i = eachindex(C) + C[i] += A[aoffs + i] * b + end end end end diff --git a/test/matmul.jl b/test/matmul.jl index d405225b..3d5a9c96 100644 --- a/test/matmul.jl +++ b/test/matmul.jl @@ -960,11 +960,23 @@ Base.:*(x::Float64, a::A32092) = x * a.x end @testset "strong zero" begin - @testset for α in Any[false, 0.0, 0], n in 1:4 - C = ones(n, n) - A = fill!(zeros(n, n), NaN) - B = ones(n, n) + @testset for α in Any[false, 0.0, 0], n in 1:4, T in (Float16, Float64) + C = ones(T, n) + A = fill(T(NaN), n, n) + B = ones(T, n) @test mul!(copy(C), A, B, α, 1.0) == C + C = ones(T, n, n) + B = ones(T, n, n) + @test mul!(copy(C), A, B, α, 1.0) == C + end + @testset for α in Any[false, 0.0, 0], β in Any[false, 0.0, 0], n in 1:4, T in (Float16, Float64) + C = fill(T(NaN), n) + A = fill(T(NaN), n, n) + B = fill(T(NaN), n) + @test iszero(mul!(copy(C), A, B, α, β)) + C = fill(T(NaN), n, n) + B = fill(T(NaN), n, n) + @test iszero(mul!(copy(C), A, B, α, β)) end end From fb7bee457e7467ef619e42dca5100795e3ea732a Mon Sep 17 00:00:00 2001 From: Lorenzo Van Munoz <66997677+lxvm@users.noreply.github.com> Date: Fri, 15 Aug 2025 05:11:02 -0700 Subject: [PATCH 08/12] Fix overflow and type issues with generic det of a triangular matrix (#1418) Fixes #1415 --------- Co-authored-by: Steven G. Johnson --- src/generic.jl | 2 +- test/generic.jl | 15 +++++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/generic.jl b/src/generic.jl index 7c69f20c..3b9a4232 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1821,7 +1821,7 @@ julia> det(BigInt[1 0; 2 2]) # exact integer determinant function det(A::AbstractMatrix{T}) where {T} if istriu(A) || istril(A) S = promote_type(T, typeof((one(T)*zero(T) + zero(T))/one(T))) - return convert(S, det(UpperTriangular(A))) + return prod(Base.Fix1(convert, S), @view A[diagind(A)]; init=one(S)) end return det(lu(A; check = false)) end diff --git a/test/generic.jl b/test/generic.jl index 11ff874f..10bdc19c 100644 --- a/test/generic.jl +++ b/test/generic.jl @@ -25,6 +25,9 @@ using .Main.FillArrays isdefined(Main, :SizedArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "SizedArrays.jl")) using .Main.SizedArrays +isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl")) +using .Main.Furlongs + Random.seed!(123) n = 5 # should be odd @@ -103,6 +106,18 @@ n = 5 # should be odd @testset "det with nonstandard Number type" begin elty <: Real && @test det(Dual.(triu(A), zero(A))) isa Dual end + if elty <: Int + @testset "det no overflow - triangular" begin + A = diagm([typemax(elty), typemax(elty)]) + @test det(A) == det(float(A)) + end + end + @testset "det with units - triangular" begin + for dim in 0:4 + A = diagm(Furlong.(ones(elty, dim))) + @test det(A) == Furlong{dim}(one(elty)) + end + end end @testset "diag" begin From 4648cc4b6e92d419f53e93ad8cd489906f17d3bb Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 18 Aug 2025 15:42:47 -0400 Subject: [PATCH 09/12] Update ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c1611fbf..f592ded7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -18,7 +18,7 @@ jobs: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 with: - version: 'nightly' + version: '1.12' - name: Generate docs run: | julia --color=yes --project=docs -e 'using Pkg; Pkg.instantiate()' From 664fa901730e8a39b068682b111477dd71f127de Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 18 Aug 2025 15:43:56 -0400 Subject: [PATCH 10/12] Update ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f592ded7..6ab3baf0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -18,7 +18,7 @@ jobs: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 with: - version: '1.12' + version: 'pre' - name: Generate docs run: | julia --color=yes --project=docs -e 'using Pkg; Pkg.instantiate()' From ebc850684e8c4f74c043fef175331d39351b82d7 Mon Sep 17 00:00:00 2001 From: Viral Shah Date: Mon, 18 Aug 2025 16:20:36 -0400 Subject: [PATCH 11/12] Use pre to run 1.12 pre-release in tests --- .buildkite/runtests.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.buildkite/runtests.yml b/.buildkite/runtests.yml index 936dc49d..ba4d5c9d 100644 --- a/.buildkite/runtests.yml +++ b/.buildkite/runtests.yml @@ -13,7 +13,7 @@ steps: gid: 1000 # Julia installation inside the sandbox - JuliaCI/julia#v1: - version: "nightly" + version: "pre" arch: "i686" command: | julia --color=yes --code-coverage=@ .ci/run_tests.jl @@ -37,7 +37,7 @@ steps: gid: 1000 # Julia installation inside the sandbox - JuliaCI/julia#v1: - version: "nightly" + version: "pre" command: | julia --color=yes --code-coverage=@ .ci/run_tests.jl agents: @@ -49,7 +49,7 @@ steps: - label: ":macos: macos-aarch64" plugins: - JuliaCI/julia#v1: - version: "nightly" + version: "pre" - JuliaCI/julia-coverage#v1: codecov: true command: | @@ -62,7 +62,7 @@ steps: - label: ":windows: windows-x86_64" plugins: - JuliaCI/julia#v1: - version: "nightly" + version: "pre" - JuliaCI/julia-coverage#v1: codecov: true command: | From e4eeff1b3217270df8d41639bbf54f41de335abb Mon Sep 17 00:00:00 2001 From: Viral Shah Date: Mon, 18 Aug 2025 16:24:44 -0400 Subject: [PATCH 12/12] Try julia version 1.12 in buildkite --- .buildkite/runtests.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.buildkite/runtests.yml b/.buildkite/runtests.yml index ba4d5c9d..bffd3c3b 100644 --- a/.buildkite/runtests.yml +++ b/.buildkite/runtests.yml @@ -13,7 +13,7 @@ steps: gid: 1000 # Julia installation inside the sandbox - JuliaCI/julia#v1: - version: "pre" + version: "1.12" arch: "i686" command: | julia --color=yes --code-coverage=@ .ci/run_tests.jl @@ -37,7 +37,7 @@ steps: gid: 1000 # Julia installation inside the sandbox - JuliaCI/julia#v1: - version: "pre" + version: "1.12" command: | julia --color=yes --code-coverage=@ .ci/run_tests.jl agents: @@ -49,7 +49,7 @@ steps: - label: ":macos: macos-aarch64" plugins: - JuliaCI/julia#v1: - version: "pre" + version: "1.12" - JuliaCI/julia-coverage#v1: codecov: true command: | @@ -62,7 +62,7 @@ steps: - label: ":windows: windows-x86_64" plugins: - JuliaCI/julia#v1: - version: "pre" + version: "1.12" - JuliaCI/julia-coverage#v1: codecov: true command: |