Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LinearAlgebra: replace some hardcoded loop ranges with axes #56243

Merged
merged 4 commits into from
Oct 25, 2024
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 46 additions & 49 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -705,18 +705,15 @@ 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
@inbounds begin
for j = 1:n
nrmj::Tsum = 0
for i = 1:m
nrmj += norm(A[i,j])
end
nrm = max(nrm,nrmj)
for j in axes(A,2)
nrmj::Tsum = 0
for i in axes(A,1)
nrmj += norm(@inbounds A[i,j])
end
nrm = max(nrm,nrmj)
end
return convert(Tnorm, nrm)
end
Expand All @@ -732,18 +729,15 @@ 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
@inbounds begin
for i = 1:m
nrmi::Tsum = 0
for j = 1:n
nrmi += norm(A[i,j])
end
nrm = max(nrm,nrmi)
for i in axes(A,1)
nrmi::Tsum = 0
for j in axes(A,2)
nrmi += norm(@inbounds A[i,j])
end
nrm = max(nrm,nrmi)
end
return convert(Tnorm, nrm)
end
Expand Down Expand Up @@ -938,7 +932,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
Expand Down Expand Up @@ -979,11 +973,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)
Expand Down Expand Up @@ -1596,10 +1590,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 = 1:n
xi, yi = x[i], y[i]
x[i] = c *xi + s*yi
y[i] = -conj(s)*xi + c*yi
for i in 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
Expand All @@ -1619,10 +1615,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 = 1:n
xi, yi = x[i], y[i]
x[i] = c *xi + s*yi
y[i] = conj(s)*xi - c*yi
for i in eachindex(x,y)
@inbounds begin
dkarrasch marked this conversation as resolved.
Show resolved Hide resolved
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
Expand All @@ -1633,18 +1631,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 in 2:n
@inbounds x[i] /= ξ1
end
ξ1/ν
end
Expand All @@ -1655,16 +1651,16 @@ 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
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 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
axpy!(-vAj, xj, Aj)
end
return A
Expand Down Expand Up @@ -1799,9 +1795,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)
Expand Down Expand Up @@ -1991,12 +1988,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 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 1:n, 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
Expand Down