diff --git a/.buildkite/runtests.yml b/.buildkite/runtests.yml index 936dc49d..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: "nightly" + 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: "nightly" + 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: "nightly" + 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: "nightly" + version: "1.12" - JuliaCI/julia-coverage#v1: codecov: true command: | 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 diff --git a/src/adjtrans.jl b/src/adjtrans.jl index d81aa3ae..ee610be7 100644 --- a/src/adjtrans.jl +++ b/src/adjtrans.jl @@ -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)) 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/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/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/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 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/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( diff --git a/test/adjtrans.jl b/test/adjtrans.jl index bdeaa697..e5654004 100644 --- a/test/adjtrans.jl +++ b/test/adjtrans.jl @@ -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 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 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 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