From b42583cfb4c5fa4e337dcec7aaaf77d019549f6b Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 19 Oct 2024 18:55:00 +0530 Subject: [PATCH 1/4] LinearAlgebra: replace some hardcoded loop ranges with axes --- stdlib/LinearAlgebra/src/generic.jl | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 6c65c49add74b..771668dfbc51f 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -710,9 +710,9 @@ function opnorm1(A::AbstractMatrix{T}) where T Tsum = promote_type(Float64, Tnorm) nrm::Tsum = 0 @inbounds begin - for j = 1:n + for j = axes(A,2) nrmj::Tsum = 0 - for i = 1:m + for i = axes(A,1) nrmj += norm(A[i,j]) end nrm = max(nrm,nrmj) @@ -737,9 +737,9 @@ function opnormInf(A::AbstractMatrix{T}) where T Tsum = promote_type(Float64, Tnorm) nrm::Tsum = 0 @inbounds begin - for i = 1:m + for i = axes(A,1) nrmi::Tsum = 0 - for j = 1:n + for j = axes(A,2) nrmi += norm(A[i,j]) end nrm = max(nrm,nrmi) @@ -1596,7 +1596,7 @@ function rotate!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - @inbounds for i = 1:n + @inbounds for i = eachindex(x) xi, yi = x[i], y[i] x[i] = c *xi + s*yi y[i] = -conj(s)*xi + c*yi @@ -1619,7 +1619,7 @@ function reflect!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - @inbounds for i = 1:n + @inbounds for i = eachindex(x,y) xi, yi = x[i], y[i] x[i] = c *xi + s*yi y[i] = conj(s)*xi - c*yi @@ -1655,13 +1655,13 @@ end Multiplies `A` in-place by a Householder reflection on the left. It is equivalent to `A .= (I - conj(τ)*[1; x[2:end]]*[1; x[2:end]]')*A`. """ @inline function reflectorApply!(x::AbstractVector, τ::Number, A::AbstractVecOrMat) - require_one_based_indexing(x) + require_one_based_indexing(x, A) m, n = size(A, 1), size(A, 2) if length(x) != m throw(DimensionMismatch(lazy"reflector has length $(length(x)), which must match the first dimension of matrix A, $m")) end m == 0 && return A - @inbounds for j = 1:n + @inbounds for j = axes(A,2) Aj, xj = view(A, 2:m, j), view(x, 2:m) vAj = conj(τ)*(A[1, j] + dot(xj, Aj)) A[1, j] -= vAj @@ -1799,9 +1799,10 @@ julia> LinearAlgebra.det_bareiss!(M) ``` """ function det_bareiss!(M) + Base.require_one_based_indexing(M) n = checksquare(M) sign, prev = Int8(1), one(eltype(M)) - for i in 1:n-1 + for i in axes(M,2)[begin:end-1] if iszero(M[i,i]) # swap with another col to make nonzero swapto = findfirst(!iszero, @view M[i,i+1:end]) isnothing(swapto) && return zero(prev) @@ -1991,12 +1992,12 @@ function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) A = Base.unalias(B, A) if uplo == 'U' LAPACK.lacpy_size_check((m1, n1), (n < m ? n : m, n)) - for j in 1:n, i in 1:min(j,m) + for j in axes(A,2), i in 1:min(j,m) @inbounds B[i,j] = A[i,j] end else # uplo == 'L' LAPACK.lacpy_size_check((m1, n1), (m, m < n ? m : n)) - for j in 1:n, i in j:m + for j in axes(A,2), i in j:m @inbounds B[i,j] = A[i,j] end end From 34b5aa0a734ab0f153a60c7237a183f2273f26ab Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 19 Oct 2024 19:00:59 +0530 Subject: [PATCH 2/4] Update copytrito! and rotate! --- stdlib/LinearAlgebra/src/generic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 771668dfbc51f..5e6a537573e0a 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -1596,7 +1596,7 @@ function rotate!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - @inbounds for i = eachindex(x) + @inbounds for i = eachindex(x,y) xi, yi = x[i], y[i] x[i] = c *xi + s*yi y[i] = -conj(s)*xi + c*yi @@ -1992,12 +1992,12 @@ function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) A = Base.unalias(B, A) if uplo == 'U' LAPACK.lacpy_size_check((m1, n1), (n < m ? n : m, n)) - for j in axes(A,2), i in 1:min(j,m) + for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] @inbounds B[i,j] = A[i,j] end else # uplo == 'L' LAPACK.lacpy_size_check((m1, n1), (m, m < n ? m : n)) - for j in axes(A,2), i in j:m + for j in axes(A,2), i in axes(A,1)[j:end] @inbounds B[i,j] = A[i,j] end end From 3752c12e1fcc22215d1dec64edbef5a3f0004387 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 19 Oct 2024 19:13:34 +0530 Subject: [PATCH 3/4] Narrower inbounds annotations in generic LinearAlgebra functions --- stdlib/LinearAlgebra/src/generic.jl | 84 ++++++++++++++--------------- 1 file changed, 41 insertions(+), 43 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 5e6a537573e0a..386636a722ade 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -709,14 +709,12 @@ function opnorm1(A::AbstractMatrix{T}) where T Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) nrm::Tsum = 0 - @inbounds begin - for j = axes(A,2) - nrmj::Tsum = 0 - for i = axes(A,1) - nrmj += norm(A[i,j]) - end - nrm = max(nrm,nrmj) + for j = axes(A,2) + nrmj::Tsum = 0 + for i = axes(A,1) + nrmj += norm(@inbounds A[i,j]) end + nrm = max(nrm,nrmj) end return convert(Tnorm, nrm) end @@ -736,14 +734,12 @@ function opnormInf(A::AbstractMatrix{T}) where T Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) nrm::Tsum = 0 - @inbounds begin - for i = axes(A,1) - nrmi::Tsum = 0 - for j = axes(A,2) - nrmi += norm(A[i,j]) - end - nrm = max(nrm,nrmi) + for i = axes(A,1) + nrmi::Tsum = 0 + for j = axes(A,2) + nrmi += norm(@inbounds A[i,j]) end + nrm = max(nrm,nrmi) end return convert(Tnorm, nrm) end @@ -938,7 +934,7 @@ function dot(x::AbstractArray, y::AbstractArray) end s = zero(dot(first(x), first(y))) for (Ix, Iy) in zip(eachindex(x), eachindex(y)) - @inbounds s += dot(x[Ix], y[Iy]) + s += dot(@inbounds(x[Ix]), @inbounds(y[Iy])) end s end @@ -979,11 +975,11 @@ function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector) s = zero(T) i₁ = first(eachindex(x)) x₁ = first(x) - @inbounds for j in eachindex(y) - yj = y[j] + for j in eachindex(y) + yj = @inbounds y[j] if !iszero(yj) - temp = zero(adjoint(A[i₁,j]) * x₁) - @simd for i in eachindex(x) + temp = zero(adjoint(@inbounds A[i₁,j]) * x₁) + @inbounds @simd for i in eachindex(x) temp += adjoint(A[i,j]) * x[i] end s += dot(temp, yj) @@ -1596,10 +1592,12 @@ function rotate!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - @inbounds for i = eachindex(x,y) - xi, yi = x[i], y[i] - x[i] = c *xi + s*yi - y[i] = -conj(s)*xi + c*yi + for i = eachindex(x,y) + @inbounds begin + xi, yi = x[i], y[i] + x[i] = c *xi + s*yi + y[i] = -conj(s)*xi + c*yi + end end return x, y end @@ -1619,10 +1617,12 @@ function reflect!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - @inbounds for i = eachindex(x,y) - xi, yi = x[i], y[i] - x[i] = c *xi + s*yi - y[i] = conj(s)*xi - c*yi + for i = eachindex(x,y) + @inbounds begin + xi, yi = x[i], y[i] + x[i] = c *xi + s*yi + y[i] = conj(s)*xi - c*yi + end end return x, y end @@ -1633,18 +1633,16 @@ end require_one_based_indexing(x) n = length(x) n == 0 && return zero(eltype(x)) - @inbounds begin - ξ1 = x[1] - normu = norm(x) - if iszero(normu) - return zero(ξ1/normu) - end - ν = T(copysign(normu, real(ξ1))) - ξ1 += ν - x[1] = -ν - for i = 2:n - x[i] /= ξ1 - end + ξ1 = @inbounds x[1] + normu = norm(x) + if iszero(normu) + return zero(ξ1/normu) + end + ν = T(copysign(normu, real(ξ1))) + ξ1 += ν + @inbounds x[1] = -ν + for i = 2:n + @inbounds x[i] /= ξ1 end ξ1/ν end @@ -1661,10 +1659,10 @@ Multiplies `A` in-place by a Householder reflection on the left. It is equivalen throw(DimensionMismatch(lazy"reflector has length $(length(x)), which must match the first dimension of matrix A, $m")) end m == 0 && return A - @inbounds for j = axes(A,2) - Aj, xj = view(A, 2:m, j), view(x, 2:m) - vAj = conj(τ)*(A[1, j] + dot(xj, Aj)) - A[1, j] -= vAj + for j = axes(A,2) + Aj, xj = @inbounds view(A, 2:m, j), view(x, 2:m) + vAj = conj(τ)*(@inbounds(A[1, j]) + dot(xj, Aj)) + @inbounds A[1, j] -= vAj axpy!(-vAj, xj, Aj) end return A From 208429dfa3f32dea875ea3dfe7742d5684d79e0a Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 24 Oct 2024 15:20:25 +0530 Subject: [PATCH 4/4] Remove unused variables --- stdlib/LinearAlgebra/src/generic.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 386636a722ade..be8558fc3fad3 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -705,13 +705,12 @@ norm(::Missing, p::Real=2) = missing # special cases of opnorm function opnorm1(A::AbstractMatrix{T}) where T require_one_based_indexing(A) - m, n = size(A) Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) nrm::Tsum = 0 - for j = axes(A,2) + for j in axes(A,2) nrmj::Tsum = 0 - for i = axes(A,1) + for i in axes(A,1) nrmj += norm(@inbounds A[i,j]) end nrm = max(nrm,nrmj) @@ -730,13 +729,12 @@ end function opnormInf(A::AbstractMatrix{T}) where T require_one_based_indexing(A) - m,n = size(A) Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) nrm::Tsum = 0 - for i = axes(A,1) + for i in axes(A,1) nrmi::Tsum = 0 - for j = axes(A,2) + for j in axes(A,2) nrmi += norm(@inbounds A[i,j]) end nrm = max(nrm,nrmi) @@ -1592,7 +1590,7 @@ function rotate!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - for i = eachindex(x,y) + for i in eachindex(x,y) @inbounds begin xi, yi = x[i], y[i] x[i] = c *xi + s*yi @@ -1617,7 +1615,7 @@ function reflect!(x::AbstractVector, y::AbstractVector, c, s) if n != length(y) throw(DimensionMismatch(lazy"x has length $(length(x)), but y has length $(length(y))")) end - for i = eachindex(x,y) + for i in eachindex(x,y) @inbounds begin xi, yi = x[i], y[i] x[i] = c *xi + s*yi @@ -1641,7 +1639,7 @@ end ν = T(copysign(normu, real(ξ1))) ξ1 += ν @inbounds x[1] = -ν - for i = 2:n + for i in 2:n @inbounds x[i] /= ξ1 end ξ1/ν @@ -1659,7 +1657,7 @@ Multiplies `A` in-place by a Householder reflection on the left. It is equivalen throw(DimensionMismatch(lazy"reflector has length $(length(x)), which must match the first dimension of matrix A, $m")) end m == 0 && return A - for j = axes(A,2) + for j in axes(A,2) Aj, xj = @inbounds view(A, 2:m, j), view(x, 2:m) vAj = conj(τ)*(@inbounds(A[1, j]) + dot(xj, Aj)) @inbounds A[1, j] -= vAj