Merge branch 'master' into sgj/evalpoly

Author: ViralBShah

.ci/create_sysimage.jl

@@ -0,0 +1,16 @@
+using PackageCompiler, Libdl
+
+function create_sysimage_linear_algebra()
+    sysimage = tempname() * "." * Libdl.dlext
+
+    if haskey(ENV, "BUILDKITE")
+        ncores = Sys.CPU_THREADS
+    else
+        ncores = ceil(Int, Sys.CPU_THREADS / 2)
+    end
+
+    withenv("JULIA_IMAGE_THREADS" => ncores) do
+        create_sysimage(["LinearAlgebra", "Test", "Distributed", "Dates", "Printf", "Random"]; sysimage_path=sysimage, incremental=false, filter_stdlibs=true)
+    end
+    return sysimage, ncores
+end

.ci/create_sysimage_and_run_docs.jl

@@ -0,0 +1,6 @@
+include("create_sysimage.jl")
+sysimage, ncores = create_sysimage_linear_algebra()
+doc_file = joinpath(@__DIR__, "..", "docs", "make.jl")
+withenv("JULIA_NUM_THREADS" => 1) do
+    run(`$(Base.julia_cmd()) --sysimage=$sysimage --project=$(dirname(doc_file)) $doc_file`)
+end

.ci/create_sysimage_and_run_tests.jl

@@ -1,17 +1,5 @@
-using PackageCompiler, Libdl
-
-sysimage = tempname() * "." * Libdl.dlext
-
-if haskey(ENV, "BUILDKITE")
-    ncores = Sys.CPU_THREADS
-else
-    ncores = ceil(Int, Sys.CPU_THREADS / 2)
-end
-
-withenv("JULIA_IMAGE_THREADS" => ncores) do
-    create_sysimage(["LinearAlgebra", "Test", "Distributed", "Dates", "Printf", "Random"]; sysimage_path=sysimage, incremental=false, filter_stdlibs=true)
-end
-
+include("create_sysimage.jl")
+sysimage, ncores = create_sysimage_linear_algebra()
 current_dir = @__DIR__
 cmd = """Base.runtests(["LinearAlgebra"]; propagate_project=true, ncores=$ncores)"""
 withenv("JULIA_NUM_THREADS" => 1) do

.github/workflows/ci.yml

@@ -21,9 +21,9 @@ jobs:
           version: 'nightly'
       - name: Generate docs
         run: |
-          julia --project --color=yes -e ''
-          julia --project --color=yes -e 'using Pkg; Pkg.respect_sysimage_versions(false); Pkg.activate("docs"); Pkg.develop(PackageSpec(path = pwd()))'
-          julia --project=docs --color=yes docs/make.jl
+          julia --color=yes --project=.ci -e 'using Pkg; Pkg.instantiate()'
+          julia --color=yes --project=docs -e 'using Pkg; Pkg.instantiate()'
+          julia --color=yes --project=.ci .ci/create_sysimage_and_run_docs.jl
         env:
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

docs/Project.toml

@@ -1,3 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[sources]
+LinearAlgebra = {path = ".."}

src/LinearAlgebra.jl

@@ -16,6 +16,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     permutedims, permuterows!, power_by_squaring, promote_rule, real, sec, sech, setindex!,
     show, similar, sin, sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc,
     typed_hcat, vec, view, zero
+import Base: AbstractArray, AbstractMatrix, Array, Matrix
 using Base: IndexLinear, promote_eltype, promote_op, print_matrix,
     @propagate_inbounds, reduce, typed_hvcat, typed_vcat, require_one_based_indexing,
     splat, BitInteger

src/adjtrans.jl

@@ -319,8 +319,8 @@ const AdjointAbsVec{T} = Adjoint{T,<:AbstractVector}
 const AdjointAbsMat{T} = Adjoint{T,<:AbstractMatrix}
 const TransposeAbsVec{T} = Transpose{T,<:AbstractVector}
 const TransposeAbsMat{T} = Transpose{T,<:AbstractMatrix}
-const AdjOrTransAbsVec{T} = AdjOrTrans{T,<:AbstractVector}
-const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix}
+const AdjOrTransAbsVec{T,V<:AbstractVector} = AdjOrTrans{T,V}
+const AdjOrTransAbsMat{T,M<:AbstractMatrix} = AdjOrTrans{T,M}
 
 # for internal use below
 wrapperop(_) = identity

src/bidiag.jl

@@ -109,8 +109,8 @@ julia> Bidiagonal(A, :L) # contains the main diagonal and first subdiagonal of A
  ⋅  ⋅  4  4
 ```
 """
-function Bidiagonal(A::AbstractMatrix, uplo::Symbol)
-    Bidiagonal(diag(A, 0), diag(A, uplo === :U ? 1 : -1), uplo)
+function (::Type{Bi})(A::AbstractMatrix, uplo::Symbol) where {Bi<:Bidiagonal}
+    Bi(diag(A, 0), diag(A, uplo === :U ? 1 : -1), uplo)
 end
 
 
@@ -161,7 +161,7 @@ end
         # we explicitly compare the possible bands as b.band may be constant-propagated
         return @inbounds A.ev[b.index]
     else
-        return diagzero(A, Tuple(_cartinds(b))...)
+        return diagzero(A, b)
     end
 end
 
@@ -220,7 +220,12 @@ promote_rule(::Type{<:Tridiagonal}, ::Type{<:Bidiagonal}) = Tridiagonal
 AbstractMatrix{T}(A::Bidiagonal) where {T} = Bidiagonal{T}(A)
 AbstractMatrix{T}(A::Bidiagonal{T}) where {T} = copy(A)
 
-convert(::Type{T}, m::AbstractMatrix) where {T<:Bidiagonal} = m isa T ? m : T(m)::T
+function convert(::Type{T}, A::AbstractMatrix) where T<:Bidiagonal
+    checksquare(A)
+    isbanded(A, -1, 1) || throw(InexactError(:convert, T, A))
+    iszero(diagview(A, 1)) ? T(A, :L) :
+        iszero(diagview(A, -1)) ? T(A, :U) : throw(InexactError(:convert, T, A))
+end
 
 similar(B::Bidiagonal, ::Type{T}) where {T} = Bidiagonal(similar(B.dv, T), similar(B.ev, T), B.uplo)
 similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = similar(B.dv, T, dims)
@@ -405,14 +410,15 @@ end
 function diag(M::Bidiagonal, n::Integer=0)
     # every branch call similar(..., ::Int) to make sure the
     # same vector type is returned independent of n
-    v = similar(M.dv, max(0, length(M.dv)-abs(n)))
+    dinds = diagind(M, n, IndexStyle(M))
+    v = similar(M.dv, length(dinds))
     if n == 0
         copyto!(v, M.dv)
     elseif (n == 1 && M.uplo == 'U') ||  (n == -1 && M.uplo == 'L')
         copyto!(v, M.ev)
     elseif -size(M,1) <= n <= size(M,1)
-        for i in eachindex(v)
-            v[i] = M[BandIndex(n,i)]
+        for i in eachindex(v, dinds)
+            @inbounds v[i] = M[BandIndex(n,i)]
         end
     end
     return v

src/blas.jl

@@ -84,7 +84,8 @@ export
     trsm!,
     trsm
 
-using ..LinearAlgebra: libblastrampoline, BlasReal, BlasComplex, BlasFloat, BlasInt, DimensionMismatch, checksquare, chkstride1
+using ..LinearAlgebra: libblastrampoline, BlasReal, BlasComplex, BlasFloat, BlasInt,
+    DimensionMismatch, checksquare, chkstride1, SingularException
 
 include("lbt.jl")
 
@@ -128,6 +129,10 @@ Set the number of threads the BLAS library should use equal to `n::Integer`.
 
 Also accepts `nothing`, in which case julia tries to guess the default number of threads.
 Passing `nothing` is discouraged and mainly exists for historical reasons.
+
+!!! note
+    Some BLAS libraries, such as Apple Accelerate, cannot be configured to use a fixed number of threads.
+    For these backends, `set_num_threads()` is a no-op. See also [`get_num_threads`](@ref).
 """
 set_num_threads(nt::Integer)::Nothing = lbt_set_num_threads(Int32(nt))
 function set_num_threads(::Nothing)
@@ -147,6 +152,10 @@ Get the number of threads the BLAS library is using.
 
 !!! compat "Julia 1.6"
     `get_num_threads` requires at least Julia 1.6.
+
+!!! note
+    Some BLAS libraries, such as Apple Accelerate, cannot be configured to use a fixed number of threads.
+    For these backends, `get_num_threads()` always returns `1`. See also [`set_num_threads`](@ref).
 """
 get_num_threads()::Int = lbt_get_num_threads()
 
@@ -1369,6 +1378,11 @@ for (fname, elty) in ((:dtrsv_,:Float64),
                 throw(DimensionMismatch(lazy"size of A is $n != length(x) = $(length(x))"))
             end
             chkstride1(A)
+            if diag == 'N'
+                for i in 1:n
+                    iszero(A[i,i]) && throw(SingularException(i))
+                end
+            end
             px, stx = vec_pointer_stride(x, ArgumentError("input vector with 0 stride is not allowed"))
             GC.@preserve x ccall((@blasfunc($fname), libblastrampoline), Cvoid,
                 (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
@@ -2217,6 +2231,11 @@ for (mmname, smname, elty) in
             end
             chkstride1(A)
             chkstride1(B)
+            if diag == 'N'
+                for i in 1:k
+                    iszero(A[i,i]) && throw(SingularException(i))
+                end
+            end
             ccall((@blasfunc($smname), libblastrampoline), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                     Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},

src/cholesky.jl

@@ -65,7 +65,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Matrix{Float64}}:
+3×3 LowerTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -305,7 +305,7 @@ function _cholpivoted!(A::AbstractMatrix, ::Type{UpperTriangular}, tol::Real, ch
     rTA = real(eltype(A))
     # checks
     Base.require_one_based_indexing(A)
-    n = LinearAlgebra.checksquare(A)
+    n = checksquare(A)
     # initialization
     piv = collect(1:n)
     dots = zeros(rTA, n)
@@ -354,7 +354,7 @@ function _cholpivoted!(A::AbstractMatrix, ::Type{LowerTriangular}, tol::Real, ch
     rTA = real(eltype(A))
     # checks
     Base.require_one_based_indexing(A)
-    n = LinearAlgebra.checksquare(A)
+    n = checksquare(A)
     # initialization
     piv = collect(1:n)
     dots = zeros(rTA, n)
@@ -530,7 +530,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Matrix{Float64}}:
+3×3 LowerTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -664,30 +664,15 @@ copy(C::CholeskyPivoted) = CholeskyPivoted(copy(C.factors), C.uplo, C.piv, C.ran
 size(C::Union{Cholesky, CholeskyPivoted}) = size(C.factors)
 size(C::Union{Cholesky, CholeskyPivoted}, d::Integer) = size(C.factors, d)
 
-function _choleskyUfactor(Cfactors, Cuplo)
-    if Cuplo === 'U'
-        return UpperTriangular(Cfactors)
-    else
-        return copy(LowerTriangular(Cfactors)')
-    end
-end
-function _choleskyLfactor(Cfactors, Cuplo)
-    if Cuplo === 'L'
-        return LowerTriangular(Cfactors)
-    else
-        return copy(UpperTriangular(Cfactors)')
-    end
-end
-
 function getproperty(C::Cholesky, d::Symbol)
     Cfactors = getfield(C, :factors)
     Cuplo    = getfield(C, :uplo)
     if d === :U
-        _choleskyUfactor(Cfactors, Cuplo)
+        UpperTriangular(Cuplo == 'U' ? Cfactors : Cfactors')
     elseif d === :L
-        _choleskyLfactor(Cfactors, Cuplo)
+        LowerTriangular(Cuplo == 'L' ? Cfactors : Cfactors')
     elseif d === :UL
-        return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
+        return (Cuplo == 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
     else
         return getfield(C, d)
     end
@@ -704,9 +689,9 @@ function getproperty(C::CholeskyPivoted{T}, d::Symbol) where {T}
     Cfactors = getfield(C, :factors)
     Cuplo    = getfield(C, :uplo)
     if d === :U
-        _choleskyUfactor(Cfactors, Cuplo)
+        UpperTriangular(Cuplo == 'U' ? Cfactors : Cfactors')
     elseif d === :L
-        _choleskyLfactor(Cfactors, Cuplo)
+        LowerTriangular(Cuplo == 'L' ? Cfactors : Cfactors')
     elseif d === :p
         return getfield(C, :piv)
     elseif d === :P
@@ -813,7 +798,7 @@ function rdiv!(B::AbstractMatrix, C::Cholesky)
     end
 end
 
-function LinearAlgebra.rdiv!(B::AbstractMatrix, C::CholeskyPivoted)
+function rdiv!(B::AbstractMatrix, C::CholeskyPivoted)
     n = size(C, 2)
     for i in 1:size(B, 1)
         permute!(view(B, i, 1:n), C.piv)
@@ -965,7 +950,7 @@ function lowrankdowndate!(C::Cholesky, v::AbstractVector)
         s = conj(v[i]/Aii)
         s2 = abs2(s)
         if s2 > 1
-            throw(LinearAlgebra.PosDefException(i))
+            throw(PosDefException(i))
         end
         c = sqrt(1 - abs2(s))