Merge branch 'master' into jishnub/blocktri_indexing

Author: dkarrasch

.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/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)

src/cholesky.jl

@@ -65,7 +65,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Adjoint{Matrix{Float64}}}:
+3×3 LowerTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -530,7 +530,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Adjoint{Matrix{Float64}}}:
+3×3 LowerTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0

src/dense.jl

@@ -682,6 +682,22 @@ cis(A::AbstractMatrix{<:Base.HWNumber}) = exp_maybe_inplace(float.(im .* A))
 exp_maybe_inplace(A::StridedMatrix{<:Union{ComplexF32, ComplexF64}}) = exp!(A)
 exp_maybe_inplace(A) = exp(A)
 
+function copytri_maybe_inplace(A::StridedMatrix, uplo, conjugate::Bool=false, diag::Bool=false)
+    copytri!(A, uplo, conjugate, diag)
+end
+function copytri_maybe_inplace(A, uplo, conjugate::Bool=false, diag::Bool=false)
+    k = Int(diag)
+    if uplo == 'U'
+        B = triu(A, 1-k)
+        triu(A, k) + (conjugate ? copy(adjoint(B)) : copy(transpose(B)))
+    elseif uplo == 'L'
+        B = tril(A, k-1)
+        tril(A, -k) + (conjugate ? copy(adjoint(B)) : copy(transpose(B)))
+    else
+        throw(ArgumentError(lazy"uplo argument must be 'U' (upper) or 'L' (lower), got $uplo"))
+    end
+end
+
 """
     ^(b::Number, A::AbstractMatrix)
 
@@ -894,7 +910,8 @@ function log(A::AbstractMatrix)
         return applydiagonal(log, A)
     elseif ishermitian(A)
         logHermA = log(Hermitian(A))
-        return ishermitian(logHermA) ? copytri!(parent(logHermA), 'U', true) : parent(logHermA)
+        PH = parent(logHermA)
+        return ishermitian(logHermA) ? copytri_maybe_inplace(PH, 'U', true) : PH
     elseif istriu(A)
         return triu!(parent(log(UpperTriangular(A))))
     elseif isreal(A)
@@ -977,7 +994,8 @@ function sqrt(A::AbstractMatrix{T}) where {T<:Union{Real,Complex}}
         return applydiagonal(sqrt, A)
     elseif ishermitian(A)
         sqrtHermA = sqrt(Hermitian(A))
-        return ishermitian(sqrtHermA) ? copytri!(parent(sqrtHermA), 'U', true) : parent(sqrtHermA)
+        PS = parent(sqrtHermA)
+        return ishermitian(sqrtHermA) ? copytri_maybe_inplace(PS, 'U', true) : PS
     elseif istriu(A)
         return triu!(parent(sqrt(UpperTriangular(A))))
     elseif isreal(A)
@@ -1100,7 +1118,8 @@ function cos(A::AbstractMatrix{<:Real})
     if isdiag(A)
         return applydiagonal(cos, A)
     elseif issymmetric(A)
-        return copytri!(parent(cos(Symmetric(A))), 'U')
+        P = parent(cos(Symmetric(A)))
+        return copytri_maybe_inplace(P, 'U')
     end
     M = im .* float.(A)
     return real(exp_maybe_inplace(M))
@@ -1109,7 +1128,8 @@ function cos(A::AbstractMatrix{<:Complex})
     if isdiag(A)
         return applydiagonal(cos, A)
     elseif ishermitian(A)
-        return copytri!(parent(cos(Hermitian(A))), 'U', true)
+        P = parent(cos(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     M = im .* float.(A)
     N = -M
@@ -1140,7 +1160,8 @@ function sin(A::AbstractMatrix{<:Real})
     if isdiag(A)
         return applydiagonal(sin, A)
     elseif issymmetric(A)
-        return copytri!(parent(sin(Symmetric(A))), 'U')
+        P = parent(sin(Symmetric(A)))
+        return copytri_maybe_inplace(P, 'U')
     end
     M = im .* float.(A)
     return imag(exp_maybe_inplace(M))
@@ -1149,7 +1170,8 @@ function sin(A::AbstractMatrix{<:Complex})
     if isdiag(A)
         return applydiagonal(sin, A)
     elseif ishermitian(A)
-        return copytri!(parent(sin(Hermitian(A))), 'U', true)
+        P = parent(sin(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     M = im .* float.(A)
     Mneg = -M
@@ -1183,8 +1205,10 @@ julia> C
 function sincos(A::AbstractMatrix{<:Real})
     if issymmetric(A)
         symsinA, symcosA = sincos(Symmetric(A))
-        sinA = copytri!(parent(symsinA), 'U')
-        cosA = copytri!(parent(symcosA), 'U')
+        Psin = parent(symsinA)
+        Pcos = parent(symcosA)
+        sinA = copytri_maybe_inplace(Psin, 'U')
+        cosA = copytri_maybe_inplace(Pcos, 'U')
         return sinA, cosA
     end
     M =  im .* float.(A)
@@ -1194,8 +1218,10 @@ end
 function sincos(A::AbstractMatrix{<:Complex})
     if ishermitian(A)
         hermsinA, hermcosA = sincos(Hermitian(A))
-        sinA = copytri!(parent(hermsinA), 'U', true)
-        cosA = copytri!(parent(hermcosA), 'U', true)
+        Psin = parent(hermsinA)
+        Pcos = parent(hermcosA)
+        sinA = copytri_maybe_inplace(Psin, 'U', true)
+        cosA = copytri_maybe_inplace(Pcos, 'U', true)
         return sinA, cosA
     end
     M = im .* float.(A)
@@ -1239,7 +1265,8 @@ function tan(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(tan, A)
     elseif ishermitian(A)
-        return copytri!(parent(tan(Hermitian(A))), 'U', true)
+        P = parent(tan(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     S, C = sincos(A)
     S /= C
@@ -1255,7 +1282,8 @@ function cosh(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(cosh, A)
     elseif ishermitian(A)
-        return copytri!(parent(cosh(Hermitian(A))), 'U', true)
+        P = parent(cosh(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     X = exp(A)
     negA = @. float(-A)
@@ -1272,7 +1300,8 @@ function sinh(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(sinh, A)
     elseif ishermitian(A)
-        return copytri!(parent(sinh(Hermitian(A))), 'U', true)
+        P = parent(sinh(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     X = exp(A)
     negA = @. float(-A)
@@ -1289,7 +1318,8 @@ function tanh(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(tanh, A)
     elseif ishermitian(A)
-        return copytri!(parent(tanh(Hermitian(A))), 'U', true)
+        P = parent(tanh(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     X = exp(A)
     negA = @. float(-A)
@@ -1332,7 +1362,8 @@ function acos(A::AbstractMatrix)
         return applydiagonal(acos, A)
     elseif ishermitian(A)
         acosHermA = acos(Hermitian(A))
-        return isa(acosHermA, Hermitian) ? copytri!(parent(acosHermA), 'U', true) : parent(acosHermA)
+        P = parent(acosHermA)
+        return isa(acosHermA, Hermitian) ? copytri_maybe_inplace(P, 'U', true) : P
     end
     SchurF = Schur{Complex}(schur(A))
     U = UpperTriangular(SchurF.T)
@@ -1365,7 +1396,8 @@ function asin(A::AbstractMatrix)
         return applydiagonal(asin, A)
     elseif ishermitian(A)
         asinHermA = asin(Hermitian(A))
-        return isa(asinHermA, Hermitian) ? copytri!(parent(asinHermA), 'U', true) : parent(asinHermA)
+        P = parent(asinHermA)
+        return isa(asinHermA, Hermitian) ? copytri_maybe_inplace(P, 'U', true) : P
     end
     SchurF = Schur{Complex}(schur(A))
     U = UpperTriangular(SchurF.T)
@@ -1397,7 +1429,8 @@ function atan(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(atan, A)
     elseif ishermitian(A)
-        return copytri!(parent(atan(Hermitian(A))), 'U', true)
+        P = parent(atan(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     SchurF = Schur{Complex}(schur(A))
     U = im * UpperTriangular(SchurF.T)
@@ -1418,7 +1451,8 @@ function acosh(A::AbstractMatrix)
         return applydiagonal(acosh, A)
     elseif ishermitian(A)
         acoshHermA = acosh(Hermitian(A))
-        return isa(acoshHermA, Hermitian) ? copytri!(parent(acoshHermA), 'U', true) : parent(acoshHermA)
+        P = parent(acoshHermA)
+        return isa(acoshHermA, Hermitian) ? copytri_maybe_inplace(P, 'U', true) : P
     end
     SchurF = Schur{Complex}(schur(A))
     U = UpperTriangular(SchurF.T)
@@ -1438,7 +1472,8 @@ function asinh(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(asinh, A)
     elseif ishermitian(A)
-        return copytri!(parent(asinh(Hermitian(A))), 'U', true)
+        P = parent(asinh(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     SchurF = Schur{Complex}(schur(A))
     U = UpperTriangular(SchurF.T)
@@ -1458,7 +1493,8 @@ function atanh(A::AbstractMatrix)
     if isdiag(A)
         return applydiagonal(atanh, A)
     elseif ishermitian(A)
-        return copytri!(parent(atanh(Hermitian(A))), 'U', true)
+        P = parent(atanh(Hermitian(A)))
+        return copytri_maybe_inplace(P, 'U', true)
     end
     SchurF = Schur{Complex}(schur(A))
     U = UpperTriangular(SchurF.T)