Merge branch 'master' into dk/mul

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

@@ -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
 

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/diagonal.jl

@@ -213,7 +213,7 @@ zeroslike(::Type{M}, sz::Tuple{Integer, Vararg{Integer}}) where {M<:AbstractMatr
     if b.band == 0
         @inbounds r = D.diag[b.index]
     else
-        r = diagzero(D, Tuple(_cartinds(b))...)
+        r = diagzero(D, b)
     end
     r
 end
@@ -353,6 +353,13 @@ function rmul!(A::AbstractMatrix, D::Diagonal)
     end
     return A
 end
+# A' = A' * D => A = D' * A
+# This uses the fact that D' is a Diagonal
+function rmul!(A::AdjOrTransAbsMat, D::Diagonal)
+    f = wrapperop(A)
+    lmul!(f(D), f(A))
+    A
+end
 # T .= T * D
 function rmul!(T::Tridiagonal, D::Diagonal)
     matmul_size_check(size(T), size(D))
@@ -388,6 +395,13 @@ function lmul!(D::Diagonal, T::Tridiagonal)
     end
     return T
 end
+# A' = D * A' => A = A * D'
+# This uses the fact that D' is a Diagonal
+function lmul!(D::Diagonal, A::AdjOrTransAbsMat)
+    f = wrapperop(A)
+    rmul!(f(A), f(D))
+    A
+end
 
 @inline function __muldiag_nonzeroalpha!(out, D::Diagonal, B, alpha::Number, beta::Number)
     @inbounds for j in axes(B, 2)

src/triangular.jl

@@ -139,10 +139,10 @@ julia> UnitUpperTriangular(A)
 """
 UnitUpperTriangular
 
-const UpperOrUnitUpperTriangular{T,S} = Union{UpperTriangular{T,S}, UnitUpperTriangular{T,S}}
-const LowerOrUnitLowerTriangular{T,S} = Union{LowerTriangular{T,S}, UnitLowerTriangular{T,S}}
-const UpperOrLowerTriangular{T,S} = Union{UpperOrUnitUpperTriangular{T,S}, LowerOrUnitLowerTriangular{T,S}}
-const UnitUpperOrUnitLowerTriangular{T,S} = Union{UnitUpperTriangular{T,S}, UnitLowerTriangular{T,S}}
+const UpperOrUnitUpperTriangular{T,S<:AbstractMatrix{T}} = Union{UpperTriangular{T,S}, UnitUpperTriangular{T,S}}
+const LowerOrUnitLowerTriangular{T,S<:AbstractMatrix{T}} = Union{LowerTriangular{T,S}, UnitLowerTriangular{T,S}}
+const UpperOrLowerTriangular{T,S<:AbstractMatrix{T}} = Union{UpperOrUnitUpperTriangular{T,S}, LowerOrUnitLowerTriangular{T,S}}
+const UnitUpperOrUnitLowerTriangular{T,S<:AbstractMatrix{T}} = Union{UnitUpperTriangular{T,S}, UnitLowerTriangular{T,S}}
 
 uppertriangular(M) = UpperTriangular(M)
 lowertriangular(M) = LowerTriangular(M)

src/tridiag.jl

@@ -474,6 +474,19 @@ end
     end
 end
 
+@inline function getindex(A::SymTridiagonal, b::BandIndex)
+    @boundscheck checkbounds(A, b)
+    if b.band == 0
+        return symmetric((@inbounds A.dv[b.index]), :U)::symmetric_type(eltype(A.dv))
+    elseif b.band == -1
+        return copy(transpose(@inbounds A.ev[b.index])) # materialized for type stability
+    elseif b.band == 1
+        return @inbounds A.ev[b.index]
+    else
+        return diagzero(A, b)
+    end
+end
+
 Base._reverse(A::SymTridiagonal, dims) = reverse!(Matrix(A); dims)
 Base._reverse(A::SymTridiagonal, dims::Colon) = SymTridiagonal(reverse(A.dv), reverse(A.ev))
 Base._reverse!(A::SymTridiagonal, dims::Colon) = (reverse!(A.dv); reverse!(A.ev); A)