Skip to content
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
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