Skip to content
Merged
Show file tree
Hide file tree
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
8 changes: 4 additions & 4 deletions .buildkite/runtests.yml
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ steps:
gid: 1000
# Julia installation inside the sandbox
- JuliaCI/julia#v1:
version: "nightly"
version: "1.12"
arch: "i686"
command: |
julia --color=yes --code-coverage=@ .ci/run_tests.jl
Expand All @@ -37,7 +37,7 @@ steps:
gid: 1000
# Julia installation inside the sandbox
- JuliaCI/julia#v1:
version: "nightly"
version: "1.12"
command: |
julia --color=yes --code-coverage=@ .ci/run_tests.jl
agents:
Expand All @@ -49,7 +49,7 @@ steps:
- label: ":macos: macos-aarch64"
plugins:
- JuliaCI/julia#v1:
version: "nightly"
version: "1.12"
- JuliaCI/julia-coverage#v1:
codecov: true
command: |
Expand All @@ -62,7 +62,7 @@ steps:
- label: ":windows: windows-x86_64"
plugins:
- JuliaCI/julia#v1:
version: "nightly"
version: "1.12"
- JuliaCI/julia-coverage#v1:
codecov: true
command: |
Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -745,7 +745,7 @@ LinearAlgebra.BLAS.trsv

```@docs
LinearAlgebra.BLAS.ger!
# xGERU
LinearAlgebra.BLAS.geru!
# xGERC
LinearAlgebra.BLAS.her!
# xHPR
Expand Down
24 changes: 24 additions & 0 deletions src/adjtrans.jl
Original file line number Diff line number Diff line change
Expand Up @@ -536,3 +536,27 @@ 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))

# 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))
2 changes: 1 addition & 1 deletion src/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
2 changes: 1 addition & 1 deletion src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 7 additions & 5 deletions src/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
22 changes: 22 additions & 0 deletions src/qr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/svd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
9 changes: 4 additions & 5 deletions src/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)"))
Expand Down Expand Up @@ -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(
Expand Down
22 changes: 22 additions & 0 deletions test/adjtrans.jl
Original file line number Diff line number Diff line change
Expand Up @@ -749,4 +749,26 @@ 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

@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
15 changes: 15 additions & 0 deletions test/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
20 changes: 16 additions & 4 deletions test/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion test/svd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down