diff --git a/docs/Project.toml b/docs/Project.toml index 8e05444..834fbd6 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,5 +4,5 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" [compat] -Documenter = "0.26" -Quaternions = "0.4" \ No newline at end of file +Documenter = "1" +Quaternions = "0.7" diff --git a/docs/src/index.md b/docs/src/index.md index 745c222..51d4e56 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,6 +7,9 @@ CurrentModule = GenericLinearAlgebra ``` ```@docs +ldlt +ldlt! +numnegevals svd! svdvals! ``` diff --git a/src/ldlt.jl b/src/ldlt.jl index a5622a8..8d54232 100644 --- a/src/ldlt.jl +++ b/src/ldlt.jl @@ -147,7 +147,7 @@ end See [`ldlt`](@ref) """ -function LinearAlgebra.ldlt!(A::Hermitian{T}, blocksize::Int = 32 ÷ sizeof(T)) where T +function LinearAlgebra.ldlt!(A::Hermitian{T}, blocksize::Int = max(1, 128 ÷ sizeof(T))) where T if A.uplo === 'U' _ldlt_upper_blocked!(A.data, blocksize) else @@ -157,7 +157,7 @@ function LinearAlgebra.ldlt!(A::Hermitian{T}, blocksize::Int = 32 ÷ sizeof(T)) end """ - ldlt(A::Hermitian)::LTLt + ldlt(A::Hermitian, blocksize::Int)::LTLt A Hermitian LDL factorization of `A` such that `A = L*D*L'` if `A.uplo == 'L'` and `A = U'*D*U` if `A.uplo == 'U'. Hence, the `t` is a bit of a misnomer, @@ -171,35 +171,37 @@ The factorization has three properties: `d`, `D`, and `L` which is respectively a vector of the diagonal elements of `D`, the `Diagonal` matrix `D` and the `L` matrix when `A.uplo == 'L'` or the `adjoint` of the `U` matrix when `A.uplo == 'U'`. -The factorization is blocked. Currently, the blocking size is set to `32 ÷ sizeof(eltype(A))` -based on very rudimentary benchmarking on my laptop. +The `blocksize` argument controls the block size in the blocked algorithm. +Currently, the blocking size is set to `128 ÷ sizeof(eltype(A))` +based on very rudimentary benchmarking on my laptop. Most users won't need +adjust this argument. # Examples ```jldoctest -julia> ldlt(Hermitian([1 1; 1 -1])) -LDLt{Int64, Hermitian{Int64, Matrix{Int64}}} +julia> ldlt(Hermitian([1//1 1; 1 -1])) +LDLt{Rational{Int64}, Hermitian{Rational{Int64}, Matrix{Rational{Int64}}}} L factor: -2×2 UnitLowerTriangular{Int64, Adjoint{Int64, Matrix{Int64}}}: +2×2 UnitLowerTriangular{Rational{Int64}, Adjoint{Rational{Int64}, Matrix{Rational{Int64}}}}: 1 ⋅ 1 1 D factor: -2×2 Diagonal{Int64, SubArray{Int64, 1, Base.ReshapedArray{Int64, 1, Hermitian{Int64, Matrix{Int64}}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}, Tuple{StepRange{Int64, Int64}}, false}}: +2×2 Diagonal{Rational{Int64}, SubArray{Rational{Int64}, 1, Base.ReshapedArray{Rational{Int64}, 1, Hermitian{Rational{Int64}, Matrix{Rational{Int64}}}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}, Tuple{StepRange{Int64, Int64}}, false}}: 1 ⋅ ⋅ -2 -julia> ldlt(Hermitian([1 1; 1 1], :L)) -LDLt{Int64, Hermitian{Int64, Matrix{Int64}}} +julia> ldlt(Hermitian([1//1 1; 1 1], :L)) +LDLt{Rational{Int64}, Hermitian{Rational{Int64}, Matrix{Rational{Int64}}}} L factor: -2×2 UnitLowerTriangular{Int64, Matrix{Int64}}: +2×2 UnitLowerTriangular{Rational{Int64}, Matrix{Rational{Int64}}}: 1 ⋅ 1 1 D factor: -2×2 Diagonal{Int64, SubArray{Int64, 1, Base.ReshapedArray{Int64, 1, Hermitian{Int64, Matrix{Int64}}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}, Tuple{StepRange{Int64, Int64}}, false}}: +2×2 Diagonal{Rational{Int64}, SubArray{Rational{Int64}, 1, Base.ReshapedArray{Rational{Int64}, 1, Hermitian{Rational{Int64}, Matrix{Rational{Int64}}}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}, Tuple{StepRange{Int64, Int64}}, false}}: 1 ⋅ ⋅ 0 ``` """ -LinearAlgebra.ldlt(A::Hermitian{T}, blocksize::Int = 32 ÷ sizeof(T)) where T = ldlt!(LinearAlgebra.copy_oftype(A, typeof(oneunit(T) / one(T)))) +LinearAlgebra.ldlt(A::Hermitian{T}, blocksize::Int = max(1, 128 ÷ sizeof(T))) where T = ldlt!(LinearAlgebra.copy_oftype(A, typeof(oneunit(T) / one(T)))) function Base.getproperty(F::LDLt{<:Any, <:Hermitian}, d::Symbol) Fdata = getfield(F, :data) diff --git a/src/svd.jl b/src/svd.jl index 013032d..07ab684 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -567,6 +567,8 @@ LinearAlgebra.svdvals!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real} = svdvals!(A [, tol]) Generic computation of singular values. + +# Examples ```jldoctest julia> using LinearAlgebra, GenericLinearAlgebra, Quaternions @@ -574,12 +576,12 @@ julia> n = 20; julia> H = [big(1)/(i + j - 1) for i in 1:n, j in 1:n]; # The Hilbert matrix -julia> Float64(svdvals(H)[end]/svdvals(Float64.(H))[end] - 1) # The relative error of the LAPACK based solution in 64 bit floating point. --0.9999999999447275 +julia> round(svdvals(H)[end]/svdvals(Float64.(H))[end] - 1, sigdigits=8) # The relative error of the LAPACK based solution rounded to eight significant digits. +-1.0 julia> A = qr([Quaternion(randn(4)...) for i in 1:3, j in 1:3]).Q * Diagonal([3, 2, 1]) * - qr([Quaternion(randn(4)...) for i in 1:3, j in 1:3]).Q'; # A quaternion matrix with the singular value 1, 2, and 3. + qr([Quaternion(randn(4)...) for i in 1:3, j in 1:3]).Q'; julia> svdvals(A) ≈ [3, 2, 1] true @@ -617,7 +619,7 @@ A generic singular value decomposition (SVD). The implementation only uses Julia ```jldoctest julia> svd(big.([1 2; 3 4])) -SVD{BigFloat, BigFloat, Matrix{BigFloat}} +SVD{BigFloat, BigFloat, Matrix{BigFloat}, Vector{BigFloat}} U factor: 2×2 Matrix{BigFloat}: -0.404554 0.914514