diff --git a/.buildkite/runtests.yml b/.buildkite/runtests.yml
index a21eb633..936dc49d 100644
--- a/.buildkite/runtests.yml
+++ b/.buildkite/runtests.yml
@@ -16,8 +16,7 @@ steps:
           version: "nightly"
           arch: "i686"
     command: |
-      julia --color=yes --project=.ci -e 'using Pkg; Pkg.instantiate()'
-      julia --color=yes --project=.ci .ci/create_sysimage_and_run_tests.jl
+      julia --color=yes --code-coverage=@ .ci/run_tests.jl
     agents:
       queue: "julia"
       os: "linux"
@@ -40,8 +39,7 @@ steps:
       - JuliaCI/julia#v1:
           version: "nightly"
     command: |
-      julia --color=yes --project=.ci --code-coverage=@ -e 'using Pkg; Pkg.instantiate()'
-      julia --color=yes --project=.ci --code-coverage=@ .ci/create_sysimage_and_run_tests.jl
+      julia --color=yes --code-coverage=@ .ci/run_tests.jl
     agents:
       queue: "julia"
       os: "linux"
@@ -55,8 +53,7 @@ steps:
       - JuliaCI/julia-coverage#v1:
           codecov: true
     command: |
-      julia --color=yes --project=.ci --code-coverage=@ -e 'using Pkg; Pkg.instantiate()'
-      julia --color=yes --project=.ci --code-coverage=@ .ci/create_sysimage_and_run_tests.jl
+      julia --color=yes --code-coverage=@ .ci/run_tests.jl
     agents:
       queue: "julia"
       os: "macos"
@@ -69,8 +66,7 @@ steps:
       - JuliaCI/julia-coverage#v1:
           codecov: true
     command: |
-      julia --color=yes --project=.ci --code-coverage=@ -e 'using Pkg; Pkg.instantiate()'
-      julia --color=yes --project=.ci --code-coverage=@ .ci/create_sysimage_and_run_tests.jl
+      julia --color=yes --code-coverage=@ .ci/run_tests.jl
     agents:
       queue: "julia"
       os: "windows"
diff --git a/.ci/Manifest.toml b/.ci/Manifest.toml
deleted file mode 100644
index 4f888b42..00000000
--- a/.ci/Manifest.toml
+++ /dev/null
@@ -1,222 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
-julia_version = "1.12.0-DEV"
-manifest_format = "2.0"
-project_hash = "b25b51369ac1c4a63619c6f77655347661ec5439"
-
-[[deps.ArgTools]]
-uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
-version = "1.1.2"
-
-[[deps.Artifacts]]
-uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
-version = "1.11.0"
-
-[[deps.Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-version = "1.11.0"
-
-[[deps.CompilerSupportLibraries_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.3.0+1"
-
-[[deps.Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-version = "1.11.0"
-
-[[deps.Distributed]]
-deps = ["Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-version = "1.11.0"
-
-[[deps.Downloads]]
-deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
-uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
-version = "1.6.0"
-
-[[deps.FileWatching]]
-uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
-version = "1.11.0"
-
-[[deps.Glob]]
-git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
-uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
-version = "1.3.1"
-
-[[deps.InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-version = "1.11.0"
-
-[[deps.JuliaSyntaxHighlighting]]
-deps = ["StyledStrings"]
-uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011"
-version = "1.12.0"
-
-[[deps.LazyArtifacts]]
-deps = ["Artifacts", "Pkg"]
-uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
-version = "1.11.0"
-
-[[deps.LibCURL]]
-deps = ["LibCURL_jll", "MozillaCACerts_jll"]
-uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.4"
-
-[[deps.LibCURL_jll]]
-deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll", "Zlib_jll", "nghttp2_jll"]
-uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "8.11.1+1"
-
-[[deps.LibGit2]]
-deps = ["LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-version = "1.11.0"
-
-[[deps.LibGit2_jll]]
-deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll"]
-uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
-version = "1.9.0+0"
-
-[[deps.LibSSH2_jll]]
-deps = ["Artifacts", "Libdl", "OpenSSL_jll"]
-uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.11.3+1"
-
-[[deps.Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-version = "1.11.0"
-
-[[deps.LinearAlgebra]]
-deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
-path = ".."
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-version = "1.12.0"
-
-[[deps.Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-version = "1.11.0"
-
-[[deps.Markdown]]
-deps = ["Base64", "JuliaSyntaxHighlighting", "StyledStrings"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-version = "1.11.0"
-
-[[deps.MozillaCACerts_jll]]
-uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2024.12.31"
-
-[[deps.NetworkOptions]]
-uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
-version = "1.3.0"
-
-[[deps.OpenBLAS_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
-uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.29+0"
-
-[[deps.OpenSSL_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.0.15+2"
-
-[[deps.PackageCompiler]]
-deps = ["Artifacts", "Glob", "LazyArtifacts", "Libdl", "Pkg", "Printf", "RelocatableFolders", "TOML", "UUIDs", "p7zip_jll"]
-git-tree-sha1 = "5d13e5b70011762b74f86fc08385303589f80272"
-uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
-version = "2.2.0"
-
-[[deps.Pkg]]
-deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.12.0"
-
-    [deps.Pkg.extensions]
-    REPLExt = "REPL"
-
-    [deps.Pkg.weakdeps]
-    REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[deps.Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-version = "1.11.0"
-
-[[deps.Random]]
-deps = ["SHA"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-version = "1.11.0"
-
-[[deps.RelocatableFolders]]
-deps = ["SHA", "Scratch"]
-git-tree-sha1 = "ffdaf70d81cf6ff22c2b6e733c900c3321cab864"
-uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
-version = "1.0.1"
-
-[[deps.SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-version = "0.7.0"
-
-[[deps.Scratch]]
-deps = ["Dates"]
-git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386"
-uuid = "6c6a2e73-6563-6170-7368-637461726353"
-version = "1.2.1"
-
-[[deps.Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-version = "1.11.0"
-
-[[deps.Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-version = "1.11.0"
-
-[[deps.StyledStrings]]
-uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
-version = "1.11.0"
-
-[[deps.TOML]]
-deps = ["Dates"]
-uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
-version = "1.0.3"
-
-[[deps.Tar]]
-deps = ["ArgTools", "SHA"]
-uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
-version = "1.10.0"
-
-[[deps.Test]]
-deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-version = "1.11.0"
-
-[[deps.UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-version = "1.11.0"
-
-[[deps.Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
-version = "1.11.0"
-
-[[deps.Zlib_jll]]
-deps = ["Libdl"]
-uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.3.1+2"
-
-[[deps.libblastrampoline_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.12.0+0"
-
-[[deps.nghttp2_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.64.0+1"
-
-[[deps.p7zip_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.5.0+2"
diff --git a/.ci/Project.toml b/.ci/Project.toml
deleted file mode 100644
index 3916fe72..00000000
--- a/.ci/Project.toml
+++ /dev/null
@@ -1,8 +0,0 @@
-[deps]
-Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
-Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
-Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/.ci/create_sysimage.jl b/.ci/create_sysimage.jl
deleted file mode 100644
index 43aea023..00000000
--- a/.ci/create_sysimage.jl
+++ /dev/null
@@ -1,16 +0,0 @@
-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
diff --git a/.ci/create_sysimage_and_run_docs.jl b/.ci/create_sysimage_and_run_docs.jl
deleted file mode 100644
index 376219bb..00000000
--- a/.ci/create_sysimage_and_run_docs.jl
+++ /dev/null
@@ -1,6 +0,0 @@
-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
diff --git a/.ci/create_sysimage_and_run_tests.jl b/.ci/create_sysimage_and_run_tests.jl
deleted file mode 100644
index 071fe1fd..00000000
--- a/.ci/create_sysimage_and_run_tests.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-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
-    run(`$(Base.julia_cmd()) --sysimage=$sysimage --project=$current_dir -e $cmd`)
-end
diff --git a/.ci/run_tests.jl b/.ci/run_tests.jl
new file mode 100644
index 00000000..879b73b1
--- /dev/null
+++ b/.ci/run_tests.jl
@@ -0,0 +1,10 @@
+if haskey(ENV, "BUILDKITE")
+    ncores = Sys.CPU_THREADS
+else
+    ncores = ceil(Int, Sys.CPU_THREADS / 2)
+end
+
+proj = abspath(joinpath(@__DIR__, ".."))
+cmd = """Base.runtests(["LinearAlgebra"]; propagate_project=true, ncores=$ncores)"""
+run(addenv(`$(Base.julia_cmd()) --project=$proj --compiled-modules=existing -e $cmd`,
+    "JULIA_NUM_THREADS" => 1, "JULIA_PRUNE_OLD_LA" => true))
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 7c423f04..c1611fbf 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -21,9 +21,8 @@ jobs:
           version: 'nightly'
       - name: Generate docs
         run: |
-          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
+          julia --color=yes --project=docs --compiled-modules=existing docs/make.jl
         env:
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
diff --git a/README.md b/README.md
index a39a7248..c0951dcc 100644
--- a/README.md
+++ b/README.md
@@ -37,10 +37,37 @@ This package performs some type piracy and is also included in the sysimage, whi
 To use a development version of this package, you can choose one of the following methods:
 
 1. **Change the UUID in the project file and load the package:**
-   This approach will produce warnings and may lead to method ambiguities between the development version and the one in the sysimage, but it can be used for basic experimentation.
+
+   Change the UUID line in `Project.toml` as
+   ```diff
+   - uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+   + uuid = "27e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+   ```
+
+   Start `julia` as
+   ```console
+   JULIA_PRUNE_OLD_LA=true julia +nightly --compiled-modules=existing --project
+   ```
+   where it is assumed that one is already within the `LinearAlgebra` directory (otherwise, adjust
+   the project path accordingly). The `julia +nightly` command above assumes that `juliaup` is being used
+   to launch `julia`, but one may substitute this with the path to the julia executable.
+
+   Within the `julia` session, load the `LinearAlgebra` module after pruning the one in the sysimage. This may be done as
+   ```julia
+   include("test/prune_old_LA.jl") && using LinearAlgebra
+   ```
+
+   Note that loading the test files in the REPL will automatically carry out the pruning to ensure that the development version of the package is used in the tests.
+
+   If you are contributing to the repo using this method, it may be convenient to ignore the local changes to `Project.toml` by running
+   ```console
+   git update-index --skip-worktree Project.toml
+   ```
 
 2. **Build Julia with the custom `LinearAlgebra` commit:**
-   Modify the commit in `stdlib/LinearAlgebra.version` and build Julia.
+
+   Modify the commit in [`stdlib/LinearAlgebra.version`](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra.version) and build Julia. This requires one to push the development branch
+   to `GitHub` or an equivalent platform.
 
 3. **Build a custom sysimage with the new `LinearAlgebra`:**
    - Install `PackageCompiler`.
diff --git a/docs/make.jl b/docs/make.jl
index 351fec43..0f8b79d9 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,3 +1,7 @@
+withenv("JULIA_PRUNE_OLD_LA" => "true") do
+    include("../test/prune_old_LA.jl")
+end
+
 using LinearAlgebra
 using Documenter: DocMeta, makedocs, deploydocs, HTML
 
diff --git a/src/LinearAlgebra.jl b/src/LinearAlgebra.jl
index 2c76dcdc..6f6e21ba 100644
--- a/src/LinearAlgebra.jl
+++ b/src/LinearAlgebra.jl
@@ -179,7 +179,8 @@ public AbstractTriangular,
         symmetric,
         symmetric_type,
         zeroslike,
-        matprod_dest
+        matprod_dest,
+        fillstored!
 
 const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
 const BlasReal = Union{Float64,Float32}
@@ -393,13 +394,13 @@ control over the factorization of `A`.
 ```jldoctest
 julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];
 
-julia> X = [1; 2.5; 3];
+julia> B = [1, 2.5, 3];
 
-julia> Y = zero(X);
+julia> Y = similar(B); # use similar since there is no need to read from it
 
-julia> ldiv!(Y, qr(A), X);
+julia> ldiv!(Y, qr(A), B); # you may also try qr!(A) to further reduce allocation
 
-julia> Y ≈ A\\X
+julia> Y ≈ A \\ B
 true
 ```
 """
@@ -425,13 +426,13 @@ control over the factorization of `A`.
 ```jldoctest
 julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];
 
-julia> X = [1; 2.5; 3];
+julia> B = [1, 2.5, 3];
 
-julia> Y = copy(X);
+julia> B0 = copy(B); # a backup copy to facilitate testing
 
-julia> ldiv!(qr(A), X);
+julia> ldiv!(lu(A), B); # you may also try lu!(A) to further reduce allocation
 
-julia> X ≈ A\\Y
+julia> B ≈ A \\ B0
 true
 ```
 """
@@ -726,6 +727,8 @@ end
 (\)(F::TransposeFactorization{T,<:LU}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} =
     ldiv(F, B)
 
+const default_peakflops_size = Int === Int32 ? 2048 : 4096
+
 """
     LinearAlgebra.peakflops(n::Integer=4096; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false)
 
@@ -750,7 +753,7 @@ of the problem that is solved on each processor.
     This function requires at least Julia 1.1. In Julia 1.0 it is available from
     the standard library `InteractiveUtils`.
 """
-function peakflops(n::Integer=4096; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false)
+function peakflops(n::Integer=default_peakflops_size; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false)
     t = zeros(Float64, ntrials)
     for i=1:ntrials
         a = ones(eltype,n,n)
diff --git a/src/bidiag.jl b/src/bidiag.jl
index 89e5711c..ca95e386 100644
--- a/src/bidiag.jl
+++ b/src/bidiag.jl
@@ -65,16 +65,13 @@ julia> Bl = Bidiagonal(dv, ev, :L) # ev is on the first subdiagonal
  ⋅  ⋅  9  4
 ```
 """
-function Bidiagonal(dv::V, ev::V, uplo::Symbol) where {T,V<:AbstractVector{T}}
-    Bidiagonal{T,V}(dv, ev, uplo)
-end
-function Bidiagonal(dv::V, ev::V, uplo::AbstractChar) where {T,V<:AbstractVector{T}}
+function Bidiagonal(dv::V, ev::V, uplo::Union{Symbol, AbstractChar}) where {T,V<:AbstractVector{T}}
     Bidiagonal{T,V}(dv, ev, uplo)
 end
 
 #To allow Bidiagonal's where the "dv" is Vector{T} and "ev" Vector{S},
 #where T and S can be promoted
-function Bidiagonal(dv::Vector{T}, ev::Vector{S}, uplo::Symbol) where {T,S}
+function Bidiagonal(dv::Vector{T}, ev::Vector{S}, uplo::Union{Symbol, AbstractChar}) where {T,S}
     TS = promote_type(T,S)
     return Bidiagonal{TS,Vector{TS}}(dv, ev, uplo)
 end
@@ -963,7 +960,7 @@ function _mul_bitrisym!(C::AbstractVecOrMat, A::Bidiagonal, B::AbstractVecOrMat,
     if A.uplo == 'U'
         u = A.ev
         @inbounds begin
-            for j = 1:nB
+            for j in axes(B,2)
                 b₀, b₊ = B[1, j], B[2, j]
                 _modify!(_add, d[1]*b₀ + u[1]*b₊, C, (1, j))
                 for i = 2:nA - 1
@@ -976,7 +973,7 @@ function _mul_bitrisym!(C::AbstractVecOrMat, A::Bidiagonal, B::AbstractVecOrMat,
     else
         l = A.ev
         @inbounds begin
-            for j = 1:nB
+            for j in axes(B,2)
                 b₀, b₊ = B[1, j], B[2, j]
                 _modify!(_add, d[1]*b₀, C, (1, j))
                 for i = 2:nA - 1
@@ -996,7 +993,7 @@ function _mul_bitrisym!(C::AbstractVecOrMat, A::TriSym, B::AbstractVecOrMat, _ad
     d = _diag(A, 0)
     u = _diag(A, 1)
     @inbounds begin
-        for j = 1:nB
+        for j in axes(B,2)
             b₀, b₊ = B[1, j], B[2, j]
             _modify!(_add, d[1]*b₀ + u[1]*b₊, C, (1, j))
             for i = 2:nA - 1
@@ -1028,7 +1025,7 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::TriSym, _add::MulAddMul)
         B21 = Bl[1]
         Bmm = Bd[m]
         Bm₋1m = Bu[m-1]
-        for i in 1:n
+        for i in axes(A,1)
             _modify!(_add, A[i,1] * B11 + A[i, 2] * B21, C, (i, 1))
             _modify!(_add, A[i, m-1] * Bm₋1m + A[i, m] * Bmm, C, (i, m))
         end
@@ -1037,7 +1034,7 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::TriSym, _add::MulAddMul)
             Bj₋1j = Bu[j-1]
             Bjj = Bd[j]
             Bj₊1j = Bl[j]
-            for i = 1:n
+            for i in axes(A,1)
                 _modify!(_add, A[i, j-1] * Bj₋1j + A[i, j]*Bjj + A[i, j+1] * Bj₊1j, C, (i, j))
             end
         end
@@ -1052,17 +1049,17 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::Bidiagonal, _add::MulAdd
     (iszero(m) || iszero(n)) && return C
     iszero(_add.alpha) && return _rmul_or_fill!(C, _add.beta)
     @inbounds if B.uplo == 'U'
-        for j in n:-1:2, i in 1:m
+        for j in reverse(axes(A,2)[2:end]), i in axes(A,1)
             _modify!(_add, A[i,j] * B.dv[j] + A[i,j-1] * B.ev[j-1], C, (i, j))
         end
-        for i in 1:m
+        for i in axes(A,1)
             _modify!(_add, A[i,1] * B.dv[1], C, (i, 1))
         end
     else # uplo == 'L'
-        for j in 1:n-1, i in 1:m
+        for j in axes(A,2)[1:end-1], i in axes(A,1)
             _modify!(_add, A[i,j] * B.dv[j] + A[i,j+1] * B.ev[j], C, (i, j))
         end
-        for i in 1:m
+        for i in axes(A,1)
             _modify!(_add, A[i,n] * B.dv[n], C, (i, n))
         end
     end
@@ -1138,7 +1135,7 @@ function _dibimul!(C::AbstractMatrix, A::Diagonal, B::Bidiagonal, _add)
             if B.uplo == 'L'
                 C[2,1] += _add(Ad[2]*Bev[1])
             end
-            for col in 2:n-1
+            for col in axes(A,1)[2:end-1]
                 evrow = col+rowshift
                 C[evrow, col] += _add(Ad[evrow]*Bev[col - evshift])
                 C[col, col] += _add(Ad[col]*Bdv[col])
@@ -1219,7 +1216,7 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
     @inbounds if B.uplo == 'U'
         x₀ = x[1]
         r = dot(x[1], dv[1], y[1])
-        for j in 2:nx-1
+        for j in eachindex(x)[2:end-1]
             x₋, x₀ = x₀, x[j]
             r += dot(adjoint(ev[j-1])*x₋ + adjoint(dv[j])*x₀, y[j])
         end
@@ -1229,7 +1226,7 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
         x₀ = x[1]
         x₊ = x[2]
         r = dot(adjoint(dv[1])*x₀ + adjoint(ev[1])*x₊, y[1])
-        for j in 2:nx-1
+        for j in eachindex(x)[2:end-1]
             x₀, x₊ = x₊, x[j+1]
             r += dot(adjoint(dv[j])*x₀ + adjoint(ev[j])*x₊, y[j])
         end
@@ -1260,15 +1257,15 @@ function ldiv!(c::AbstractVecOrMat, A::Bidiagonal, b::AbstractVecOrMat)
     zi = findfirst(iszero, A.dv)
     isnothing(zi) || throw(SingularException(zi))
 
-    @inbounds for j in 1:nb
+    @inbounds for j in axes(b,2)
         if A.uplo == 'L' #do colwise forward substitution
             c[1,j] = bi1 = A.dv[1] \ b[1,j]
-            for i in 2:N
+            for i in eachindex(A.dv)[2:end]
                 c[i,j] = bi1 = A.dv[i] \ (b[i,j] - A.ev[i - 1] * bi1)
             end
         else #do colwise backward substitution
             c[N,j] = bi1 = A.dv[N] \ b[N,j]
-            for i in (N - 1):-1:1
+            for i in reverse(eachindex(A.ev))
                 c[i,j] = bi1 = A.dv[i] \ (b[i,j] - A.ev[i] * bi1)
             end
         end
diff --git a/src/cholesky.jl b/src/cholesky.jl
index 03f7c273..2952503e 100644
--- a/src/cholesky.jl
+++ b/src/cholesky.jl
@@ -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)
@@ -813,7 +813,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 +965,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))
 
diff --git a/src/dense.jl b/src/dense.jl
index 234391e5..234b1022 100644
--- a/src/dense.jl
+++ b/src/dense.jl
@@ -297,6 +297,9 @@ Return a view into the `k`th diagonal of the matrix `M`.
 
 See also [`diag`](@ref), [`diagind`](@ref).
 
+!!! compat "Julia 1.12"
+    This function requires Julia 1.12 or later.
+
 # Examples
 ```jldoctest
 julia> A = [1 2 3; 4 5 6; 7 8 9]
@@ -399,6 +402,7 @@ Construct a matrix with elements of the vector as diagonal elements.
 By default, the matrix is square and its size is given by
 `length(v)`, but a non-square size `m`×`n` can be specified
 by passing `m,n` as the first arguments.
+The diagonal will be zero-padded if necessary.
 
 # Examples
 ```jldoctest
@@ -407,6 +411,13 @@ julia> diagm([1,2,3])
  1  0  0
  0  2  0
  0  0  3
+
+julia> diagm(4, 5, [1,2,3])
+4×5 Matrix{Int64}:
+ 1  0  0  0  0
+ 0  2  0  0  0
+ 0  0  3  0  0
+ 0  0  0  0  0
 ```
 """
 diagm(v::AbstractVector) = diagm(0 => v)
@@ -890,8 +901,12 @@ julia> log(A)
 """
 function log(A::AbstractMatrix)
     # If possible, use diagonalization
-    if isdiag(A)
-        return applydiagonal(log, A)
+    if isdiag(A) && eltype(A) <: Union{Real,Complex}
+        if eltype(A) <: Real && all(>=(0), diagview(A))
+            return applydiagonal(log, A)
+        else
+            return applydiagonal(log∘complex, A)
+        end
     elseif ishermitian(A)
         logHermA = log(Hermitian(A))
         return ishermitian(logHermA) ? copytri!(parent(logHermA), 'U', true) : parent(logHermA)
diff --git a/src/diagonal.jl b/src/diagonal.jl
index c4f88276..ddee8981 100644
--- a/src/diagonal.jl
+++ b/src/diagonal.jl
@@ -743,16 +743,16 @@ end
 kron(A::Diagonal, B::Diagonal) = Diagonal(kron(A.diag, B.diag))
 
 function kron(A::Diagonal, B::SymTridiagonal)
-    kdv = kron(diag(A), B.dv)
+    kdv = kron(A.diag, B.dv)
     # We don't need to drop the last element
-    kev = kron(diag(A), _pushzero(_evview(B)))
+    kev = kron(A.diag, _pushzero(_evview(B)))
     SymTridiagonal(kdv, kev)
 end
 function kron(A::Diagonal, B::Tridiagonal)
     # `_droplast!` is only guaranteed to work with `Vector`
-    kd = convert(Vector, kron(diag(A), B.d))
-    kdl = _droplast!(convert(Vector, kron(diag(A), _pushzero(B.dl))))
-    kdu = _droplast!(convert(Vector, kron(diag(A), _pushzero(B.du))))
+    kd = convert(Vector, kron(A.diag, B.d))
+    kdl = _droplast!(convert(Vector, kron(A.diag, _pushzero(B.dl))))
+    kdu = _droplast!(convert(Vector, kron(A.diag, _pushzero(B.du))))
     Tridiagonal(kdl, kd, kdu)
 end
 
diff --git a/src/generic.jl b/src/generic.jl
index 2b03b249..7c69f20c 100644
--- a/src/generic.jl
+++ b/src/generic.jl
@@ -446,7 +446,8 @@ julia> triu(a,-3)
  1.0  1.0  1.0  1.0
 ```
 """
-function triu(M::AbstractMatrix, k::Integer = 0)
+triu(M::AbstractMatrix, k::Integer = 0) = _triu(M, Val(haszero(eltype(M))), k)
+function _triu(M::AbstractMatrix, ::Val{true}, k::Integer)
     d = similar(M)
     A = triu!(d,k)
     if iszero(k)
@@ -459,6 +460,14 @@ function triu(M::AbstractMatrix, k::Integer = 0)
     end
     return A
 end
+function _triu(M::AbstractMatrix, ::Val{false}, k::Integer)
+    d = similar(M)
+    # since the zero would need to be evaluated from the elements,
+    # we copy the array to avoid undefined references in triu!
+    copy!(d, M)
+    A = triu!(d,k)
+    return A
+end
 
 """
     tril(M, k::Integer = 0)
@@ -489,7 +498,8 @@ julia> tril(a,-3)
  1.0  0.0  0.0  0.0
 ```
 """
-function tril(M::AbstractMatrix,k::Integer=0)
+tril(M::AbstractMatrix,k::Integer=0) = _tril(M, Val(haszero(eltype(M))), k)
+function _tril(M::AbstractMatrix, ::Val{true}, k::Integer)
     d = similar(M)
     A = tril!(d,k)
     if iszero(k)
@@ -502,6 +512,14 @@ function tril(M::AbstractMatrix,k::Integer=0)
     end
     return A
 end
+function _tril(M::AbstractMatrix, ::Val{false}, k::Integer)
+    d = similar(M)
+    # since the zero would need to be evaluated from the elements,
+    # we copy the array to avoid undefined references in tril!
+    copy!(d, M)
+    A = tril!(d,k)
+    return A
+end
 
 """
     triu!(M)
@@ -1116,6 +1134,8 @@ Matrix inverse. Computes matrix `N` such that
 Computed by solving the left-division
 `N = M \\ I`.
 
+A [`SingularException`](@ref) is thrown if `M` fails numerical inversion.
+
 # Examples
 ```jldoctest
 julia> M = [2 5; 1 3]
@@ -1486,6 +1506,24 @@ function _isbanded_impl(A, kl, ku)
     beyond ku, where the elements should all be zero. The reason we separate this from the
     third group is that we may loop over all the rows using A[:, col] instead of A[rowrange, col],
     which is usually faster.
+
+    E.g., in the following 6x10 matrix with (kl,ku) = (-1,1):
+     1  1  0  0  0  0  0  0  0  0
+     1  2  2  0  0  0  0  0  0  0
+     0  2  3  3  0  0  0  0  0  0
+     0  0  3  4  4  0  0  0  0  0
+     0  0  0  4  5  5  0  0  0  0
+     0  0  0  0  5  6  6  0  0  0
+
+    last_col_nonzeroblocks: 7, as every column beyond this is entirely zero
+    last_col_emptytoprows: 2, as there are zeros above the stored bands beyond this column
+    last_col_nonemptybottomrows: 4, as there are no zeros below the stored bands beyond this column
+    colrange_onlybottomrows: 1:2, as these columns only have zeros below the stored bands
+    colrange_topbottomrows: 3:4, as these columns have zeros both above and below the stored bands
+    colrange_onlytoprows_nonzero: 5:7, as these columns only have zeros above the stored bands
+    colrange_zero_block: 8:10, as every column in this range is filled with zeros
+
+    These are used to determine which rows to check for zeros in each column.
     =#
 
     last_col_nonzeroblocks = size(A,1) + ku # fully zero rectangular block beyond this column
@@ -1493,7 +1531,9 @@ function _isbanded_impl(A, kl, ku)
     last_col_nonemptybottomrows = size(A,1) + kl - 1 # empty bottom rows after this column
 
     colrange_onlybottomrows = firstindex(A,2):min(last_col_nonemptybottomrows, last_col_emptytoprows)
-    colrange_topbottomrows = max(last_col_emptytoprows, last(colrange_onlybottomrows))+1:last_col_nonzeroblocks
+    col_topbotrows_start = max(last_col_emptytoprows, last(colrange_onlybottomrows))+1
+    col_topbotrows_end = min(last_col_nonemptybottomrows, last_col_nonzeroblocks)
+    colrange_topbottomrows = col_topbotrows_start:col_topbotrows_end
     colrange_onlytoprows_nonzero = last(colrange_topbottomrows)+1:last_col_nonzeroblocks
     colrange_zero_block = last_col_nonzeroblocks+1:lastindex(A,2)
 
diff --git a/src/hessenberg.jl b/src/hessenberg.jl
index 056b97e3..158e81c8 100644
--- a/src/hessenberg.jl
+++ b/src/hessenberg.jl
@@ -103,8 +103,8 @@ Base._reverse(A::UpperHessenberg, dims) = reverse!(Matrix(A); dims)
 
 Base.@propagate_inbounds function setindex!(A::UpperHessenberg, x, i::Integer, j::Integer)
     if i > j+1
-        x == 0 || throw(ArgumentError("cannot set index in the lower triangular part " *
-            lazy"($i, $j) of an UpperHessenberg matrix to a nonzero value ($x)"))
+        iszero(x) || throw(ArgumentError(LazyString("cannot set index in the lower triangular part ",
+            lazy"($i, $j) of an UpperHessenberg matrix to a nonzero value ($x)")))
     else
         A.data[i,j] = x
     end
diff --git a/src/lapack.jl b/src/lapack.jl
index a0c1c785..20af7ca8 100644
--- a/src/lapack.jl
+++ b/src/lapack.jl
@@ -6130,10 +6130,11 @@ for (ormhr, elty) in
             chkside(side)
             chktrans(trans)
             n = checksquare(A)
+            ntau = length(tau)
             mC, nC = size(C, 1), size(C, 2)
 
-            if n - length(tau) != 1
-                throw(DimensionMismatch(lazy"tau has length $(length(tau)), needs $(n - 1)"))
+            if max(n - 1, 0) != ntau
+                throw(DimensionMismatch(lazy"tau has length $ntau, needs $(n - 1)"))
             end
             if (side == 'L' && mC != n) || (side == 'R' && nC != n)
                 throw(DimensionMismatch("A and C matrices are not conformable"))
diff --git a/src/matmul.jl b/src/matmul.jl
index 84a83d3a..10a4c81c 100644
--- a/src/matmul.jl
+++ b/src/matmul.jl
@@ -412,52 +412,60 @@ lmul!(A, B)
 
 _vec_or_mat_str(s::Tuple{Any}) = "vector"
 _vec_or_mat_str(s::Tuple{Any,Any}) = "matrix"
-@noinline function matmul_size_check(sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}})
+function matmul_size_check(sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}})
     szA2 = get(sizeA, 2, 1)
     if szA2 != sizeB[1]
-        strA = _vec_or_mat_str(sizeA)
-        strB = _vec_or_mat_str(sizeB)
-        B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB
-        size_or_len_str_B = B_size_len isa Integer ? "length" : "size"
-        dim_or_len_str_B = B_size_len isa Integer ? "length" : "first dimension"
-        pos_str_A = LazyString(length(sizeA) == length(sizeB) ? "first " : "", strA)
-        pos_str_B = LazyString(length(sizeA) == length(sizeB) ? "second " : "", strB)
-        throw(DimensionMismatch(
-            LazyString(
-                lazy"incompatible dimensions for matrix multiplication: ",
-                lazy"tried to multiply a $strA of size $sizeA with a $strB of $size_or_len_str_B $B_size_len. ",
-                lazy"The second dimension of the $pos_str_A: $szA2, does not match the $dim_or_len_str_B of the $pos_str_B: $(sizeB[1])."
-            )
-            )
-        )
+        matmul_size_check_error(sizeA, sizeB)
     end
     return nothing
 end
-@noinline function matmul_size_check(sizeC::Tuple{Integer,Vararg{Integer}}, sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}})
+@noinline function matmul_size_check_error(sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}})
+    strA = _vec_or_mat_str(sizeA)
+    strB = _vec_or_mat_str(sizeB)
+    szA2 = get(sizeA, 2, 1)
+    B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB
+    size_or_len_str_B = B_size_len isa Integer ? "length" : "size"
+    dim_or_len_str_B = B_size_len isa Integer ? "length" : "first dimension"
+    pos_str_A = LazyString(length(sizeA) == length(sizeB) ? "first " : "", strA)
+    pos_str_B = LazyString(length(sizeA) == length(sizeB) ? "second " : "", strB)
+    throw(DimensionMismatch(
+        LazyString(
+            "incompatible dimensions for matrix multiplication: ",
+            lazy"tried to multiply a $strA of size $sizeA with a $strB of $size_or_len_str_B $B_size_len. ",
+            lazy"The second dimension of the $pos_str_A: $szA2, does not match the $dim_or_len_str_B of the $pos_str_B: $(sizeB[1])."
+        )
+        )
+    )
+end
+function matmul_size_check(sizeC::Tuple{Integer,Vararg{Integer}}, sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}})
     matmul_size_check(sizeA, sizeB)
     szB2 = get(sizeB, 2, 1)
     szC2 = get(sizeC, 2, 1)
     if sizeC[1] != sizeA[1] || szC2 != szB2
-        strA = _vec_or_mat_str(sizeA)
-        strB = _vec_or_mat_str(sizeB)
-        strC = _vec_or_mat_str(sizeC)
-        C_size_len = length(sizeC) == 1 ? sizeC[1] : sizeC
-        size_or_len_str_C = C_size_len isa Integer ? "length" : "size"
-        B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB
-        size_or_len_str_B = B_size_len isa Integer ? "length" : "size"
-        destsize = length(sizeB) == length(sizeC) == 1 ? sizeA[1] : (sizeA[1], szB2)
-        size_or_len_str_dest = destsize isa Integer ? "length" : "size"
-        throw(DimensionMismatch(
-                LazyString(
-                "incompatible destination size: ",
-                lazy"the destination $strC of $size_or_len_str_C $C_size_len is incomatible with the product of a $strA of size $sizeA and a $strB of $size_or_len_str_B $B_size_len. ",
-                lazy"The destination must be of $size_or_len_str_dest $destsize."
-               )
-            )
-        )
+        matmul_size_check_error(sizeC, sizeA, sizeB)
     end
     return nothing
 end
+@noinline function matmul_size_check_error(sizeC::Tuple{Integer,Vararg{Integer}}, sizeA::Tuple{Integer,Vararg{Integer}}, sizeB::Tuple{Integer,Vararg{Integer}})
+    strA = _vec_or_mat_str(sizeA)
+    strB = _vec_or_mat_str(sizeB)
+    strC = _vec_or_mat_str(sizeC)
+    szB2 = get(sizeB, 2, 1)
+    C_size_len = length(sizeC) == 1 ? sizeC[1] : sizeC
+    size_or_len_str_C = C_size_len isa Integer ? "length" : "size"
+    B_size_len = length(sizeB) == 1 ? sizeB[1] : sizeB
+    size_or_len_str_B = B_size_len isa Integer ? "length" : "size"
+    destsize = length(sizeB) == length(sizeC) == 1 ? sizeA[1] : (sizeA[1], szB2)
+    size_or_len_str_dest = destsize isa Integer ? "length" : "size"
+    throw(DimensionMismatch(
+            LazyString(
+            "incompatible destination size: ",
+            lazy"the destination $strC of $size_or_len_str_C $C_size_len is incomatible with the product of a $strA of size $sizeA and a $strB of $size_or_len_str_B $B_size_len. ",
+            lazy"The destination must be of $size_or_len_str_dest $destsize."
+            )
+        )
+    )
+end
 
 # We may inline the matmul2x2! and matmul3x3! calls for `α == true`
 # to simplify the @stable_muladdmul branches
@@ -692,7 +700,7 @@ end
 # the aggressive constprop pushes tA and tB into gemm_wrapper!, which is needed for wrap calls within it
 # to be concretely inferred
 Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T},
-        alpha::Number, beta::Number) where {T<:BlasFloat}
+        α::Number, β::Number) where {T<:BlasFloat}
     nC = checksquare(C)
     tA_uc = uppercase(tA) # potentially convert a WrapperChar to a Char
     if tA_uc == 'T'
@@ -708,8 +716,8 @@ Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::Abst
 
     # BLAS.syrk! only updates symmetric C
     # alternatively, make non-zero β a show-stopper for BLAS.syrk!
-    if iszero(beta) || issymmetric(C)
-        α, β = promote(alpha, beta, zero(T))
+    if iszero(β) || issymmetric(C)
+        alpha, beta = promote(α, β, zero(T))
         if (alpha isa Union{Bool,T} &&
                 beta isa Union{Bool,T} &&
                 stride(A, 1) == stride(C, 1) == 1 &&
@@ -717,7 +725,7 @@ Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::Abst
             return copytri!(BLAS.syrk!('U', tA, alpha, A, beta, C), 'U')
         end
     end
-    return gemm_wrapper!(C, tA, tAt, A, A, alpha, beta)
+    return gemm_wrapper!(C, tA, tAt, A, A, α, β)
 end
 # legacy method
 syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T}, _add::MulAddMul = MulAddMul()) where {T<:BlasFloat} =
@@ -743,7 +751,7 @@ Base.@constprop :aggressive function herk_wrapper!(C::Union{StridedMatrix{T}, St
     # Result array does not need to be initialized as long as beta==0
     #    C = Matrix{T}(undef, mA, mA)
 
-    if iszero(β) || issymmetric(C)
+    if iszero(β) || ishermitian(C)
         alpha, beta = promote(α, β, zero(T))
         if (alpha isa Union{Bool,T} &&
                 beta isa Union{Bool,T} &&
@@ -962,7 +970,7 @@ function __generic_matvecmul!(::typeof(identity), C::AbstractVector, A::Abstract
         end
         for k = eachindex(B)
             aoffs = (k-1)*Astride
-            b = @stable_muladdmul MulAddMul(alpha,beta)(B[k])
+            b = @stable_muladdmul MulAddMul(alpha,false)(B[k])
             for i = eachindex(C)
                 C[i] += A[aoffs + i] * b
             end
diff --git a/src/special.jl b/src/special.jl
index 58863858..1327af67 100644
--- a/src/special.jl
+++ b/src/special.jl
@@ -320,10 +320,33 @@ _diag_or_value(A::Diagonal) = A.diag
 _diag_or_value(A::UniformScaling) = A.λ
 
 # fill[stored]! methods
+"""
+    fillstored!(A::AbstractMatrix, x)
+
+Fill only the stored indices of a structured matrix `A` with the value `x`.
+
+# Example
+```jldoctest
+julia> A = Tridiagonal(zeros(2), zeros(3), zeros(2))
+3×3 Tridiagonal{Float64, Vector{Float64}}:
+ 0.0  0.0   ⋅
+ 0.0  0.0  0.0
+  ⋅   0.0  0.0
+
+julia> LinearAlgebra.fillstored!(A, 2)
+3×3 Tridiagonal{Float64, Vector{Float64}}:
+ 2.0  2.0   ⋅
+ 2.0  2.0  2.0
+  ⋅   2.0  2.0
+```
+"""
 fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A)
 fillstored!(A::Bidiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
 fillstored!(A::Tridiagonal, x) = (fill!(A.dl, x); fill!(A.d, x); fill!(A.du, x); A)
-fillstored!(A::SymTridiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
+function fillstored!(A::SymTridiagonal, x)
+    issymmetric(x) || throw(ArgumentError("cannot set a diagonal entry of a SymTridiagonal to an asymmetric value"))
+    (fill!(A.dv, x); fill!(A.ev, x); A)
+end
 
 _small_enough(A::Union{Diagonal, Bidiagonal}) = size(A, 1) <= 1
 _small_enough(A::Tridiagonal) = size(A, 1) <= 2
diff --git a/src/structuredbroadcast.jl b/src/structuredbroadcast.jl
index 2c1e7bfa..cef4a788 100644
--- a/src/structuredbroadcast.jl
+++ b/src/structuredbroadcast.jl
@@ -84,9 +84,9 @@ function structured_broadcast_alloc(bc, ::Type{Bidiagonal}, ::Type{ElType}, n) w
     return Bidiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n1), uplo)
 end
 structured_broadcast_alloc(bc, ::Type{SymTridiagonal}, ::Type{ElType}, n) where {ElType} =
-    SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n-1))
+    SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, max(0,n-1)))
 structured_broadcast_alloc(bc, ::Type{Tridiagonal}, ::Type{ElType}, n) where {ElType} =
-    Tridiagonal(Array{ElType}(undef, n-1),Array{ElType}(undef, n),Array{ElType}(undef, n-1))
+    Tridiagonal(Array{ElType}(undef, max(0,n-1)),Array{ElType}(undef, n),Array{ElType}(undef, max(0,n-1)))
 structured_broadcast_alloc(bc, ::Type{LowerTriangular}, ::Type{ElType}, n) where {ElType} =
     LowerTriangular(Array{ElType}(undef, n, n))
 structured_broadcast_alloc(bc, ::Type{UpperTriangular}, ::Type{ElType}, n) where {ElType} =
@@ -264,13 +264,24 @@ function copyto!(dest::Tridiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
     return dest
 end
 
+# Recursively replace wrapped matrices by their parents to improve broadcasting performance
+# We may do this because the indexing within `copyto!` is restricted to the stored indices
+preprocess_broadcasted(::Type{T}, A) where {T} = _preprocess_broadcasted(T, A)
+function preprocess_broadcasted(::Type{T}, bc::Broadcasted) where {T}
+    args = map(x -> preprocess_broadcasted(T, x), bc.args)
+    Broadcast.Broadcasted(bc.f, args, bc.axes)
+end
+_preprocess_broadcasted(::Type{LowerTriangular}, A) = lowertridata(A)
+_preprocess_broadcasted(::Type{UpperTriangular}, A) = uppertridata(A)
+
 function copyto!(dest::LowerTriangular, bc::Broadcasted{<:StructuredMatrixStyle})
     isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc))
     axs = axes(dest)
     axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
+    bc_unwrapped = preprocess_broadcasted(LowerTriangular, bc)
     for j in axs[2]
         for i in j:axs[1][end]
-            @inbounds dest.data[i,j] = bc[CartesianIndex(i, j)]
+            @inbounds dest.data[i,j] = bc_unwrapped[CartesianIndex(i, j)]
         end
     end
     return dest
@@ -280,9 +291,10 @@ function copyto!(dest::UpperTriangular, bc::Broadcasted{<:StructuredMatrixStyle}
     isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc))
     axs = axes(dest)
     axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
+    bc_unwrapped = preprocess_broadcasted(UpperTriangular, bc)
     for j in axs[2]
         for i in 1:j
-            @inbounds dest.data[i,j] = bc[CartesianIndex(i, j)]
+            @inbounds dest.data[i,j] = bc_unwrapped[CartesianIndex(i, j)]
         end
     end
     return dest
diff --git a/src/symmetric.jl b/src/symmetric.jl
index 11b34a99..1edb0570 100644
--- a/src/symmetric.jl
+++ b/src/symmetric.jl
@@ -401,13 +401,11 @@ fill!(A::HermOrSym, x) = fillstored!(A, x)
 function fillstored!(A::HermOrSym{T}, x) where T
     xT = convert(T, x)
     if isa(A, Hermitian)
-        isreal(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a nonreal value"))
-    end
-    if A.uplo == 'U'
-        fillband!(A.data, xT, 0, size(A,2)-1)
-    else # A.uplo == 'L'
-        fillband!(A.data, xT, 1-size(A,1), 0)
+        ishermitian(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a non-hermitian value"))
+    elseif isa(A, Symmetric)
+        issymmetric(xT) || throw(ArgumentError("cannot fill Symmetric matrix with an asymmetric value"))
     end
+    applytri(A -> fillstored!(A, xT), A)
     return A
 end
 
diff --git a/src/triangular.jl b/src/triangular.jl
index 4a03cff6..6ce33dc7 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -599,14 +599,15 @@ end
     copytrito!(dest, U, U isa UpperOrUnitUpperTriangular ? 'L' : 'U')
     return dest
 end
-@propagate_inbounds function _copyto!(dest::StridedMatrix, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix})
+@propagate_inbounds function _copyto!(dest::StridedMatrix,
+    U::Union{UnitUpperOrUnitLowerTriangular, UpperOrLowerTriangular{<:Any, <:StridedMatrix}})
     U2 = Base.unalias(dest, U)
     copyto_unaliased!(dest, U2)
     return dest
 end
 # for strided matrices, we explicitly loop over the arrays to improve cache locality
 # This fuses the copytrito! for the two halves
-@inline function copyto_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular{<:Any, <:StridedMatrix})
+@inline function copyto_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular)
     @boundscheck checkbounds(dest, axes(U)...)
     isunit = U isa UnitUpperTriangular
     for col in axes(dest,2)
@@ -619,7 +620,7 @@ end
     end
     return dest
 end
-@inline function copyto_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular{<:Any, <:StridedMatrix})
+@inline function copyto_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular)
     @boundscheck checkbounds(dest, axes(L)...)
     isunit = L isa UnitLowerTriangular
     for col in axes(dest,2)
@@ -658,10 +659,8 @@ end
 
 uppertridata(A) = A
 lowertridata(A) = A
-# we restrict these specializations only to strided matrices to avoid cases where an UpperTriangular type
-# doesn't share its indexing with the parent
-uppertridata(A::UpperTriangular{<:Any, <:StridedMatrix}) = parent(A)
-lowertridata(A::LowerTriangular{<:Any, <:StridedMatrix}) = parent(A)
+uppertridata(A::UpperTriangular) = parent(A)
+lowertridata(A::LowerTriangular) = parent(A)
 
 @inline _rscale_add!(A::AbstractTriangular, B::AbstractTriangular, C::Number, alpha::Number, beta::Number) =
     @stable_muladdmul _triscale!(A, B, C, MulAddMul(alpha, beta))
@@ -1272,6 +1271,15 @@ end
 # Generic routines #
 ####################
 
+function _set_diag!(B::UpperOrLowerTriangular, x)
+    # get a mutable array to modify the diagonal
+    Bm = parent(B) isa StridedArray ? B : copy!(similar(B), B)
+    for i in diagind(Bm.data, IndexStyle(Bm.data))
+        Bm.data[i] = x
+    end
+    Bm
+end
+
 for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
                    (LowerTriangular, UnitLowerTriangular))
     tstrided = t{<:Any, <:StridedMaybeAdjOrTransMat}
@@ -1285,10 +1293,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (*)(A::$unitt, x::Number)
             B = $t(A.data)*x
-            for i in axes(A, 1)
-                B.data[i,i] = x
-            end
-            return B
+            _set_diag!(B, oneunit(eltype(A)) * x)
         end
 
         (*)(x::Number, A::$t) = $t(x*A.data)
@@ -1300,10 +1305,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (*)(x::Number, A::$unitt)
             B = x*$t(A.data)
-            for i in axes(A, 1)
-                B.data[i,i] = x
-            end
-            return B
+            _set_diag!(B, x * oneunit(eltype(A)))
         end
 
         (/)(A::$t, x::Number) = $t(A.data/x)
@@ -1315,11 +1317,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (/)(A::$unitt, x::Number)
             B = $t(A.data)/x
-            invx = inv(x)
-            for i in axes(A, 1)
-                B.data[i,i] = invx
-            end
-            return B
+            _set_diag!(B, oneunit(eltype(A)) / x)
         end
 
         (\)(x::Number, A::$t) = $t(x\A.data)
@@ -1331,11 +1329,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (\)(x::Number, A::$unitt)
             B = x\$t(A.data)
-            invx = inv(x)
-            for i in axes(A, 1)
-                B.data[i,i] = invx
-            end
-            return B
+            _set_diag!(B, x \ oneunit(eltype(A)))
         end
     end
 end
diff --git a/test/abstractq.jl b/test/abstractq.jl
index 302f2abf..ff1499ae 100644
--- a/test/abstractq.jl
+++ b/test/abstractq.jl
@@ -2,6 +2,8 @@
 
 module TestAbstractQ
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test
 using LinearAlgebra
 using LinearAlgebra: AbstractQ, AdjointQ
diff --git a/test/addmul.jl b/test/addmul.jl
index 903e3b17..f6b49d4e 100644
--- a/test/addmul.jl
+++ b/test/addmul.jl
@@ -2,6 +2,8 @@
 
 module TestAddmul
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Base: rtoldefault
 using Test
 using LinearAlgebra
diff --git a/test/adjtrans.jl b/test/adjtrans.jl
index 50565583..bdeaa697 100644
--- a/test/adjtrans.jl
+++ b/test/adjtrans.jl
@@ -2,6 +2,8 @@
 
 module TestAdjointTranspose
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
diff --git a/test/ambiguous_exec.jl b/test/ambiguous_exec.jl
index 7b89c0a4..09693791 100644
--- a/test/ambiguous_exec.jl
+++ b/test/ambiguous_exec.jl
@@ -1,5 +1,9 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
+module TestAmbiguity
+
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra
 let ambig = detect_ambiguities(LinearAlgebra; recursive=true)
     @test isempty(ambig)
@@ -19,3 +23,5 @@ let ambig = detect_ambiguities(LinearAlgebra; recursive=true)
     @test isempty(expect)
     @test good
 end
+
+end # module
diff --git a/test/bidiag.jl b/test/bidiag.jl
index a39fa027..982eaf2d 100644
--- a/test/bidiag.jl
+++ b/test/bidiag.jl
@@ -2,6 +2,8 @@
 
 module TestBidiagonal
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasReal, BlasFloat
 
@@ -506,6 +508,10 @@ Random.seed!(1)
     @test Matrix{ComplexF64}(BD) == BD
 end
 
+@testset "Constructors with Char uplo" begin
+    @test Bidiagonal(Int8[1,2], [1], 'U') == Bidiagonal(Int8[1,2], [1], :U)
+end
+
 # Issue 10742 and similar
 let A = Bidiagonal([1,2,3], [0,0], :U)
     @test istril(A)
diff --git a/test/bitarray.jl b/test/bitarray.jl
index 76fc21a6..10b9d4b8 100644
--- a/test/bitarray.jl
+++ b/test/bitarray.jl
@@ -1,3 +1,9 @@
+# This file is a part of Julia. License is MIT: https://julialang.org/license
+
+module TestBitArray
+
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using LinearAlgebra, Test, Random
 
 tc(r1::NTuple{N,Any}, r2::NTuple{N,Any}) where {N} = all(x->tc(x...), [zip(r1,r2)...])
@@ -93,3 +99,5 @@ b2 = bitrand(v1)
 
 b1 = bitrand(n1, n1)
 @check_bit_operation diag(b1)
+
+end # module
diff --git a/test/blas.jl b/test/blas.jl
index ff1d8bab..b7f4a03a 100644
--- a/test/blas.jl
+++ b/test/blas.jl
@@ -2,6 +2,8 @@
 
 module TestBLAS
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasReal, BlasComplex
 using Libdl: dlsym, dlopen
diff --git a/test/bunchkaufman.jl b/test/bunchkaufman.jl
index 68c519d1..cf4fa488 100644
--- a/test/bunchkaufman.jl
+++ b/test/bunchkaufman.jl
@@ -2,6 +2,8 @@
 
 module TestBunchKaufman
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted
 using Base: getproperty
diff --git a/test/cholesky.jl b/test/cholesky.jl
index 6ba72432..4c1cb14b 100644
--- a/test/cholesky.jl
+++ b/test/cholesky.jl
@@ -2,6 +2,8 @@
 
 module TestCholesky
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted,
     PosDefException, RankDeficientException, chkfullrank
diff --git a/test/dense.jl b/test/dense.jl
index 1d4b330c..d0cfa007 100644
--- a/test/dense.jl
+++ b/test/dense.jl
@@ -2,6 +2,8 @@
 
 module TestDense
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal
 using Test: GenericArray
@@ -1372,4 +1374,25 @@ end
     @test !any(LinearAlgebra.getstructure(A))
 end
 
+@testset "triu/tril for block matrices" begin
+    O = ones(2,2)
+    Z = zero(O)
+    M = fill(O, 3, 3)
+    res = fill(Z, size(M))
+    res[1,2] = res[1,3] = res[2,3] = O
+    @test triu(GenericArray(M),1) == res
+    @test tril(GenericArray(M),-1) == res'
+end
+
+@testset "log for diagonal" begin
+    D = diagm([-2.0, 2.0])
+    @test log(D) ≈ log(UpperTriangular(D))
+    D = diagm([-2.0, 0.0])
+    @test log(D) ≈ log(UpperTriangular(D))
+    D = diagm([2.0, 2.0])
+    @test log(D) ≈ log(UpperTriangular(D))
+    D = diagm([2.0, 2.0*im])
+    @test log(D) ≈ log(UpperTriangular(D))
+end
+
 end # module TestDense
diff --git a/test/diagonal.jl b/test/diagonal.jl
index d1f46f47..f83ad392 100644
--- a/test/diagonal.jl
+++ b/test/diagonal.jl
@@ -2,6 +2,8 @@
 
 module TestDiagonal
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasFloat, BlasComplex
 
diff --git a/test/eigen.jl b/test/eigen.jl
index a82c7454..67e58d34 100644
--- a/test/eigen.jl
+++ b/test/eigen.jl
@@ -2,6 +2,8 @@
 
 module TestEigen
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, UtiAUi!
 
diff --git a/test/factorization.jl b/test/factorization.jl
index f80c5197..b215c246 100644
--- a/test/factorization.jl
+++ b/test/factorization.jl
@@ -1,6 +1,9 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
 module TestFactorization
+
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra
 
 @testset "equality for factorizations - $f" for f in Any[
diff --git a/test/generic.jl b/test/generic.jl
index 6d11ec82..11ff874f 100644
--- a/test/generic.jl
+++ b/test/generic.jl
@@ -2,6 +2,8 @@
 
 module TestGeneric
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using Test: GenericArray
 using LinearAlgebra: isbanded
@@ -514,17 +516,17 @@ end
 
 @testset "generic functions for checking whether matrices have banded structure" begin
     pentadiag = [1 2 3; 4 5 6; 7 8 9]
-    tridiag   = [1 2 0; 4 5 6; 0 8 9]
-    tridiagG  = GenericArray([1 2 0; 4 5 6; 0 8 9])
+    tridiag   = diagm(-1=>1:6, 1=>1:6)
+    tridiagG  = GenericArray(tridiag)
     Tridiag   = Tridiagonal(tridiag)
     ubidiag   = [1 2 0; 0 5 6; 0 0 9]
-    ubidiagG  = GenericArray([1 2 0; 0 5 6; 0 0 9])
+    ubidiagG  = GenericArray(ubidiag)
     uBidiag   = Bidiagonal(ubidiag, :U)
     lbidiag   = [1 0 0; 4 5 0; 0 8 9]
-    lbidiagG  = GenericArray([1 0 0; 4 5 0; 0 8 9])
+    lbidiagG  = GenericArray(lbidiag)
     lBidiag   = Bidiagonal(lbidiag, :L)
     adiag     = [1 0 0; 0 5 0; 0 0 9]
-    adiagG    = GenericArray([1 0 0; 0 5 0; 0 0 9])
+    adiagG    = GenericArray(adiag)
     aDiag     = Diagonal(adiag)
     @testset "istriu" begin
         @test !istriu(pentadiag)
@@ -618,6 +620,17 @@ end
             end
         end
     end
+
+    tridiag   = diagm(-1=>1:6, 1=>1:6)
+    A = [tridiag zeros(size(tridiag,1), 2)]
+    G = GenericArray(A)
+    @testset for (kl,ku) in Iterators.product(-10:10, -10:10)
+        @test isbanded(A, kl, ku) == isbanded(G, kl, ku)
+    end
+    @testset for k in -10:10
+        @test istriu(A,k) == istriu(G,k)
+        @test istril(A,k) == istril(G,k)
+    end
 end
 
 @testset "missing values" begin
diff --git a/test/givens.jl b/test/givens.jl
index 3726afa7..fb09ec4e 100644
--- a/test/givens.jl
+++ b/test/givens.jl
@@ -2,6 +2,8 @@
 
 module TestGivens
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: Givens, Rotation, givensAlgorithm
 
diff --git a/test/hessenberg.jl b/test/hessenberg.jl
index ebab3f29..8485c014 100644
--- a/test/hessenberg.jl
+++ b/test/hessenberg.jl
@@ -2,6 +2,8 @@
 
 module TestHessenberg
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
@@ -283,5 +285,12 @@ end
         @test Q2'Q2 ≈ I
     end
 end
-    
+
+@testset "multiplication with empty HessenbergQ" begin
+    @test ones(2, 0)*hessenberg(zeros(0,0)).Q == zeros(2,0)
+    @test_throws DimensionMismatch ones(2, 1)*hessenberg(zeros(0,0)).Q
+    @test hessenberg(zeros(0,0)).Q * ones(0, 2) == zeros(0,2)
+    @test_throws DimensionMismatch hessenberg(zeros(0,0)).Q * ones(1, 2)
+end
+
 end # module TestHessenberg
diff --git a/test/lapack.jl b/test/lapack.jl
index b4e3afe3..fdf5b6bd 100644
--- a/test/lapack.jl
+++ b/test/lapack.jl
@@ -2,6 +2,8 @@
 
 module TestLAPACK
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasInt
 
@@ -677,7 +679,7 @@ end
     @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
         local n = 10
         a = randn(elty, n, n)
-        A = a'*a
+        A = a'*a + I
         B = rand(elty, n, n)
         D = copy(A)
         C = copy(B)
diff --git a/test/ldlt.jl b/test/ldlt.jl
index 51abf310..3454ef43 100644
--- a/test/ldlt.jl
+++ b/test/ldlt.jl
@@ -2,6 +2,8 @@
 
 module TestLDLT
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 Random.seed!(123)
diff --git a/test/lq.jl b/test/lq.jl
index af7ab600..0d55f84f 100644
--- a/test/lq.jl
+++ b/test/lq.jl
@@ -2,6 +2,8 @@
 
 module TestLQ
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!
 
diff --git a/test/lu.jl b/test/lu.jl
index f5b4e247..ad4faa16 100644
--- a/test/lu.jl
+++ b/test/lu.jl
@@ -2,6 +2,8 @@
 
 module TestLU
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: ldiv!, BlasReal, BlasInt, BlasFloat, rdiv!
 
diff --git a/test/matmul.jl b/test/matmul.jl
index 86c75ae5..d405225b 100644
--- a/test/matmul.jl
+++ b/test/matmul.jl
@@ -2,6 +2,8 @@
 
 module TestMatmul
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Base: rtoldefault
 using Test, LinearAlgebra, Random
 using LinearAlgebra: mul!, Symmetric, Hermitian
diff --git a/test/pinv.jl b/test/pinv.jl
index c7268865..ecee397a 100644
--- a/test/pinv.jl
+++ b/test/pinv.jl
@@ -2,6 +2,8 @@
 
 module TestPinv
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 Random.seed!(12345)
diff --git a/test/prune_old_LA.jl b/test/prune_old_LA.jl
new file mode 100644
index 00000000..7e64cfa3
--- /dev/null
+++ b/test/prune_old_LA.jl
@@ -0,0 +1,79 @@
+let
+    methods_to_delete =
+    [
+    :adjoint
+    :transpose
+    :inv
+    :literal_pow
+    :\
+    :/
+    :isapprox
+    :copyto!
+    :*
+    :muladd
+    :copyto!
+    :isone
+    :kron!
+    :kron
+    :^
+    :exp
+    :cis
+    :log
+    :sqrt
+    :cbrt
+    :inv
+    :cos
+    :sin
+    :sincos
+    :tan
+    :cosh
+    :sinh
+    :tanh
+    :acos
+    :asin
+    :atan
+    :acosh
+    :asinh
+    :atanh
+    :sec
+    :sech
+    :csc
+    :csch
+    :cot
+    :coth
+    :asec
+    :asech
+    :acsc
+    :acot
+    :acoth
+    :acsch
+    ]
+
+    prune_old_LA = parse(Bool, get(ENV, "JULIA_PRUNE_OLD_LA", "false"))
+    LinalgSysImg = Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra")
+    LA = get(Base.loaded_modules, LinalgSysImg, nothing)
+    if LA !== nothing && prune_old_LA
+        @assert hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}})
+        for methss in methods_to_delete
+            meths = getglobal(Base, methss)
+            for meth in methods(meths)
+                if meth.module === LA
+                    Base.delete_method(meth)
+                end
+            end
+        end
+    end
+end
+
+# check in a separate block to ensure that the latest world age is used
+let
+    prune_old_LA = parse(Bool, get(ENV, "JULIA_PRUNE_OLD_LA", "false"))
+    LinalgSysImg = Base.PkgId(Base.UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"), "LinearAlgebra")
+    LA = get(Base.loaded_modules, LinalgSysImg, nothing)
+    if LA !== nothing && prune_old_LA
+        @assert !hasmethod(*, Tuple{Matrix{Float64}, Matrix{Float64}})
+    end
+    prune_old_LA && Base.unreference_module(LinalgSysImg)
+end
+
+pruned_old_LA = true
diff --git a/test/qr.jl b/test/qr.jl
index b6e9ce3a..2b121ca4 100644
--- a/test/qr.jl
+++ b/test/qr.jl
@@ -2,6 +2,8 @@
 
 module TestQR
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, rmul!, lmul!
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a64884dd..0466b047 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,7 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
+
+include("prune_old_LA.jl")
+
 using Test, LinearAlgebra
 
 for file in readlines(joinpath(@__DIR__, "testgroups"))
diff --git a/test/schur.jl b/test/schur.jl
index f3d494fb..e0a16b80 100644
--- a/test/schur.jl
+++ b/test/schur.jl
@@ -2,6 +2,8 @@
 
 module TestSchur
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted
 
diff --git a/test/special.jl b/test/special.jl
index ac76279f..0cb85d56 100644
--- a/test/special.jl
+++ b/test/special.jl
@@ -2,6 +2,8 @@
 
 module TestSpecial
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: rmul!, BandIndex
 
@@ -833,4 +835,28 @@ end
     end
 end
 
+@testset "fillstored!" begin
+    dv, ev = zeros(4), zeros(3)
+    D = Diagonal(dv)
+    LinearAlgebra.fillstored!(D, 2)
+    @test D == diagm(fill(2, length(dv)))
+
+    dv .= 0
+    B = Bidiagonal(dv, ev, :U)
+    LinearAlgebra.fillstored!(B, 2)
+    @test B == diagm(0=>fill(2, length(dv)), 1=>fill(2, length(ev)))
+
+    dv .= 0
+    ev .= 0
+    T = Tridiagonal(ev, dv, ev)
+    LinearAlgebra.fillstored!(T, 2)
+    @test T == diagm(-1=>fill(2, length(ev)), 0=>fill(2, length(dv)), 1=>fill(2, length(ev)))
+
+    dv .= 0
+    ev .= 0
+    ST = SymTridiagonal(dv, ev)
+    LinearAlgebra.fillstored!(ST, 2)
+    @test ST == diagm(-1=>fill(2, length(ev)), 0=>fill(2, length(dv)), 1=>fill(2, length(ev)))
+end
+
 end # module TestSpecial
diff --git a/test/structuredbroadcast.jl b/test/structuredbroadcast.jl
index 5b5cc0cd..6521ed58 100644
--- a/test/structuredbroadcast.jl
+++ b/test/structuredbroadcast.jl
@@ -1,6 +1,9 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
 module TestStructuredBroadcast
+
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
@@ -8,160 +11,164 @@ isdefined(Main, :SizedArrays) || @eval Main include(joinpath($(BASE_TEST_PATH),
 using .Main.SizedArrays
 
 @testset "broadcast[!] over combinations of scalars, structured matrices, and dense vectors/matrices" begin
-    N = 10
-    s = rand()
-    fV = rand(N)
-    fA = rand(N, N)
-    Z = copy(fA)
-    D = Diagonal(rand(N))
-    B = Bidiagonal(rand(N), rand(N - 1), :U)
-    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
-    S = SymTridiagonal(rand(N), rand(N - 1))
-    U = UpperTriangular(rand(N,N))
-    L = LowerTriangular(rand(N,N))
-    M = Matrix(rand(N,N))
-    structuredarrays = (D, B, T, U, L, M, S)
-    fstructuredarrays = map(Array, structuredarrays)
-    for (X, fX) in zip(structuredarrays, fstructuredarrays)
-        @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX))
-        @test broadcast!(sin, Z, X) == broadcast(sin, fX)
-        @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX))
-        @test broadcast!(cos, Z, X) == broadcast(cos, fX)
-        @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX))
-        @test broadcast!(*, Z, s, X) == broadcast(*, s, fX)
-        @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX))
-        @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX)
-        @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX))
-        @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX)
-
-        @test X .* 2.0 == X .* (2.0,) == fX .* 2.0
-        @test X .* 2.0 isa typeof(X)
-        @test X .* (2.0,) isa typeof(X)
-        @test isequal(X .* Inf, fX .* Inf)
-
-        two = 2
-        @test X .^ 2 ==  X .^ (2,) == fX .^ 2 == X .^ two
-        @test X .^ 2 isa typeof(X)
-        @test X .^ (2,) isa typeof(X)
-        @test X .^ two isa typeof(X)
-        @test X .^ 0 == fX .^ 0
-        @test X .^ -1 == fX .^ -1
-
-        for (Y, fY) in zip(structuredarrays, fstructuredarrays)
-            @test broadcast(+, X, Y) == broadcast(+, fX, fY)
-            @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
-            @test broadcast(*, X, Y) == broadcast(*, fX, fY)
-            @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
+    @testset for N in (0,1,2,10) # some edge cases, and a structured case
+        s = rand()
+        fV = rand(N)
+        fA = rand(N, N)
+        Z = copy(fA)
+        D = Diagonal(rand(N))
+        B = Bidiagonal(rand(N), rand(max(0,N-1)), :U)
+        T = Tridiagonal(rand(max(0,N-1)), rand(N), rand(max(0,N-1)))
+        S = SymTridiagonal(rand(N), rand(max(0,N-1)))
+        U = UpperTriangular(rand(N,N))
+        L = LowerTriangular(rand(N,N))
+        M = Matrix(rand(N,N))
+        structuredarrays = (D, B, T, U, L, M, S)
+        fstructuredarrays = map(Array, structuredarrays)
+        @testset "$(nameof(typeof(X)))" for (X, fX) in zip(structuredarrays, fstructuredarrays)
+            @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX))
+            @test broadcast!(sin, Z, X) == broadcast(sin, fX)
+            @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX))
+            @test broadcast!(cos, Z, X) == broadcast(cos, fX)
+            @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX))
+            @test broadcast!(*, Z, s, X) == broadcast(*, s, fX)
+            @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX))
+            @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX)
+            @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX))
+            @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX)
+
+            @test X .* 2.0 == X .* (2.0,) == fX .* 2.0
+            @test X .* 2.0 isa typeof(X)
+            @test X .* (2.0,) isa typeof(X)
+            @test isequal(X .* Inf, fX .* Inf)
+
+            two = 2
+            @test X .^ 2 ==  X .^ (2,) == fX .^ 2 == X .^ two
+            @test X .^ 2 isa typeof(X)
+            @test X .^ (2,) isa typeof(X)
+            @test X .^ two isa typeof(X)
+            @test X .^ 0 == fX .^ 0
+            @test X .^ -1 == fX .^ -1
+
+            for (Y, fY) in zip(structuredarrays, fstructuredarrays)
+                @test broadcast(+, X, Y) == broadcast(+, fX, fY)
+                @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
+                @test broadcast(*, X, Y) == broadcast(*, fX, fY)
+                @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
+            end
         end
-    end
-    diagonals = (D, B, T)
-    fdiagonals = map(Array, diagonals)
-    for (X, fX) in zip(diagonals, fdiagonals)
-        for (Y, fY) in zip(diagonals, fdiagonals)
-            @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY)
-            @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
-            @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY)
-            @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
+        diagonals = (D, B, T)
+        fdiagonals = map(Array, diagonals)
+        for (X, fX) in zip(diagonals, fdiagonals)
+            for (Y, fY) in zip(diagonals, fdiagonals)
+                @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY)
+                @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
+                @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY)
+                @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
+            end
         end
-    end
-    UU = UnitUpperTriangular(rand(N,N))
-    UL = UnitLowerTriangular(rand(N,N))
-    unittriangulars = (UU, UL)
-    Ttris = typeof.((UpperTriangular(parent(UU)), LowerTriangular(parent(UU))))
-    funittriangulars = map(Array, unittriangulars)
-    for (X, fX, Ttri) in zip(unittriangulars, funittriangulars, Ttris)
-        @test (Q = broadcast(sin, X); typeof(Q) == Ttri && Q == broadcast(sin, fX))
-        @test broadcast!(sin, Z, X) == broadcast(sin, fX)
-        @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX))
-        @test broadcast!(cos, Z, X) == broadcast(cos, fX)
-        @test (Q = broadcast(*, s, X); typeof(Q) == Ttri && Q == broadcast(*, s, fX))
-        @test broadcast!(*, Z, s, X) == broadcast(*, s, fX)
-        @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX))
-        @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX)
-        @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX))
-        @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX)
-
-        @test X .* 2.0 == X .* (2.0,) == fX .* 2.0
-        @test X .* 2.0 isa Ttri
-        @test X .* (2.0,) isa Ttri
-        @test isequal(X .* Inf, fX .* Inf)
-
-        two = 2
-        @test X .^ 2 ==  X .^ (2,) == fX .^ 2 == X .^ two
-        @test X .^ 2 isa typeof(X) # special cased, as isstructurepreserving
-        @test X .^ (2,) isa Ttri
-        @test X .^ two isa Ttri
-        @test X .^ 0 == fX .^ 0
-        @test X .^ -1 == fX .^ -1
-
-        for (Y, fY) in zip(unittriangulars, funittriangulars)
-            @test broadcast(+, X, Y) == broadcast(+, fX, fY)
-            @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
-            @test broadcast(*, X, Y) == broadcast(*, fX, fY)
-            @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
+        UU = UnitUpperTriangular(rand(N,N))
+        UL = UnitLowerTriangular(rand(N,N))
+        unittriangulars = (UU, UL)
+        Ttris = typeof.((UpperTriangular(parent(UU)), LowerTriangular(parent(UU))))
+        funittriangulars = map(Array, unittriangulars)
+        for (X, fX, Ttri) in zip(unittriangulars, funittriangulars, Ttris)
+            @test (Q = broadcast(sin, X); typeof(Q) == Ttri && Q == broadcast(sin, fX))
+            @test broadcast!(sin, Z, X) == broadcast(sin, fX)
+            @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX))
+            @test broadcast!(cos, Z, X) == broadcast(cos, fX)
+            @test (Q = broadcast(*, s, X); typeof(Q) == Ttri && Q == broadcast(*, s, fX))
+            @test broadcast!(*, Z, s, X) == broadcast(*, s, fX)
+            @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX))
+            @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX)
+            @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX))
+            @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX)
+
+            @test X .* 2.0 == X .* (2.0,) == fX .* 2.0
+            @test X .* 2.0 isa Ttri
+            @test X .* (2.0,) isa Ttri
+            @test isequal(X .* Inf, fX .* Inf)
+
+            two = 2
+            @test X .^ 2 ==  X .^ (2,) == fX .^ 2 == X .^ two
+            @test X .^ 2 isa typeof(X) # special cased, as isstructurepreserving
+            @test X .^ (2,) isa Ttri
+            @test X .^ two isa Ttri
+            @test X .^ 0 == fX .^ 0
+            @test X .^ -1 == fX .^ -1
+
+            for (Y, fY) in zip(unittriangulars, funittriangulars)
+                @test broadcast(+, X, Y) == broadcast(+, fX, fY)
+                @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
+                @test broadcast(*, X, Y) == broadcast(*, fX, fY)
+                @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
+            end
         end
-    end
 
-    @testset "type-stability in Bidiagonal" begin
-        B2 = @inferred (B -> .- B)(B)
-        @test B2 isa Bidiagonal
-        @test B2 == -1 * B
-        B2 = @inferred (B -> B .* 2)(B)
-        @test B2 isa Bidiagonal
-        @test B2 == B + B
-        B2 = @inferred (B -> 2 .* B)(B)
-        @test B2 isa Bidiagonal
-        @test B2 == B + B
-        B2 = @inferred (B -> B ./ 1)(B)
-        @test B2 isa Bidiagonal
-        @test B2 == B
-        B2 = @inferred (B -> 1 .\ B)(B)
-        @test B2 isa Bidiagonal
-        @test B2 == B
+        @testset "type-stability in Bidiagonal" begin
+            B2 = @inferred (B -> .- B)(B)
+            @test B2 isa Bidiagonal
+            @test B2 == -1 * B
+            B2 = @inferred (B -> B .* 2)(B)
+            @test B2 isa Bidiagonal
+            @test B2 == B + B
+            B2 = @inferred (B -> 2 .* B)(B)
+            @test B2 isa Bidiagonal
+            @test B2 == B + B
+            B2 = @inferred (B -> B ./ 1)(B)
+            @test B2 isa Bidiagonal
+            @test B2 == B
+            B2 = @inferred (B -> 1 .\ B)(B)
+            @test B2 isa Bidiagonal
+            @test B2 == B
+        end
     end
 end
 
 @testset "broadcast! where the destination is a structured matrix" begin
-    N = 5
-    A = rand(N, N)
-    sA = A + copy(A')
-    D = Diagonal(rand(N))
-    Bu = Bidiagonal(rand(N), rand(N - 1), :U)
-    Bl = Bidiagonal(rand(N), rand(N - 1), :L)
-    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
-    ◣ = LowerTriangular(rand(N,N))
-    ◥ = UpperTriangular(rand(N,N))
-    M = Matrix(rand(N,N))
-
-    @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D))
-    @test broadcast!(sin, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U)
-    @test broadcast!(sin, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L)
-    @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T))
-    @test broadcast!(sin, copy(◣), ◣) == LowerTriangular(sin.(◣))
-    @test broadcast!(sin, copy(◥), ◥) == UpperTriangular(sin.(◥))
-    @test broadcast!(sin, copy(M), M) == Matrix(sin.(M))
-    @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A))
-    @test broadcast!(*, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U)
-    @test broadcast!(*, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L)
-    @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A))
-    @test broadcast!(*, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A))
-    @test broadcast!(*, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A))
-    @test broadcast!(*, copy(M), M, A) == Matrix(broadcast(*, M, A))
-
-    @test_throws ArgumentError broadcast!(cos, copy(D), D) == Diagonal(sin.(D))
-    @test_throws ArgumentError broadcast!(cos, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U)
-    @test_throws ArgumentError broadcast!(cos, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L)
-    @test_throws ArgumentError broadcast!(cos, copy(T), T) == Tridiagonal(sin.(T))
-    @test_throws ArgumentError broadcast!(cos, copy(◣), ◣) == LowerTriangular(sin.(◣))
-    @test_throws ArgumentError broadcast!(cos, copy(◥), ◥) == UpperTriangular(sin.(◥))
-    @test_throws ArgumentError broadcast!(+, copy(D), D, A) == Diagonal(broadcast(*, D, A))
-    @test_throws ArgumentError broadcast!(+, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U)
-    @test_throws ArgumentError broadcast!(+, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L)
-    @test_throws ArgumentError broadcast!(+, copy(T), T, A) == Tridiagonal(broadcast(*, T, A))
-    @test_throws ArgumentError broadcast!(+, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A))
-    @test_throws ArgumentError broadcast!(+, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A))
-    @test_throws ArgumentError broadcast!(*, copy(◥), ◣, 2)
-    @test_throws ArgumentError broadcast!(*, copy(Bu), Bl, 2)
+    @testset for N in (0,1,2,5)
+        A = rand(N, N)
+        sA = A + copy(A')
+        D = Diagonal(rand(N))
+        Bu = Bidiagonal(rand(N), rand(max(0,N-1)), :U)
+        Bl = Bidiagonal(rand(N), rand(max(0,N-1)), :L)
+        T = Tridiagonal(rand(max(0,N-1)), rand(N), rand(max(0,N-1)))
+        ◣ = LowerTriangular(rand(N,N))
+        ◥ = UpperTriangular(rand(N,N))
+        M = Matrix(rand(N,N))
+
+        @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D))
+        @test broadcast!(sin, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U)
+        @test broadcast!(sin, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L)
+        @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T))
+        @test broadcast!(sin, copy(◣), ◣) == LowerTriangular(sin.(◣))
+        @test broadcast!(sin, copy(◥), ◥) == UpperTriangular(sin.(◥))
+        @test broadcast!(sin, copy(M), M) == Matrix(sin.(M))
+        @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A))
+        @test broadcast!(*, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U)
+        @test broadcast!(*, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L)
+        @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A))
+        @test broadcast!(*, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A))
+        @test broadcast!(*, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A))
+        @test broadcast!(*, copy(M), M, A) == Matrix(broadcast(*, M, A))
+
+        if N > 2
+            @test_throws ArgumentError broadcast!(cos, copy(D), D)
+            @test_throws ArgumentError broadcast!(cos, copy(Bu), Bu)
+            @test_throws ArgumentError broadcast!(cos, copy(Bl), Bl)
+            @test_throws ArgumentError broadcast!(cos, copy(T), T)
+            @test_throws ArgumentError broadcast!(cos, copy(◣), ◣)
+            @test_throws ArgumentError broadcast!(cos, copy(◥), ◥)
+            @test_throws ArgumentError broadcast!(+, copy(D), D, A)
+            @test_throws ArgumentError broadcast!(+, copy(Bu), Bu, A)
+            @test_throws ArgumentError broadcast!(+, copy(Bl), Bl, A)
+            @test_throws ArgumentError broadcast!(+, copy(T), T, A)
+            @test_throws ArgumentError broadcast!(+, copy(◣), ◣, A)
+            @test_throws ArgumentError broadcast!(+, copy(◥), ◥, A)
+            @test_throws ArgumentError broadcast!(*, copy(◥), ◣, 2)
+            @test_throws ArgumentError broadcast!(*, copy(Bu), Bl, 2)
+        end
+    end
 end
 
 @testset "map[!] over combinations of structured matrices" begin
@@ -381,4 +388,12 @@ end
     @test ind == CartesianIndex(1,1)
 end
 
+@testset "nested triangular broadcast" begin
+    for T in (LowerTriangular, UpperTriangular)
+        L = T(rand(Int,4,4))
+        M = Matrix(L)
+        @test L .+ L .+ 0 .+ L .+ 0 .- L == 2M
+    end
+end
+
 end
diff --git a/test/svd.jl b/test/svd.jl
index 33c02517..fdcfcb30 100644
--- a/test/svd.jl
+++ b/test/svd.jl
@@ -2,6 +2,8 @@
 
 module TestSVD
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted
 
diff --git a/test/symmetric.jl b/test/symmetric.jl
index 68f0ff47..17500537 100644
--- a/test/symmetric.jl
+++ b/test/symmetric.jl
@@ -2,6 +2,8 @@
 
 module TestSymmetric
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
@@ -1163,4 +1165,13 @@ end
     end
 end
 
+@testset "fillstored!" begin
+    A = zeros(4,4)
+    for T in (Symmetric, Hermitian)
+        A .= 0
+        LinearAlgebra.fillstored!(T(A), 2)
+        @test all(==(2), T(A))
+    end
+end
+
 end # module TestSymmetric
diff --git a/test/symmetriceigen.jl b/test/symmetriceigen.jl
index 7a2b88a5..36b69d29 100644
--- a/test/symmetriceigen.jl
+++ b/test/symmetriceigen.jl
@@ -2,6 +2,8 @@
 
 module TestSymmetricEigen
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra
 using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
 
diff --git a/test/triangular.jl b/test/triangular.jl
index eb6814b9..1a713747 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -2,8 +2,11 @@
 
 module TestTriangular
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: errorbounds, transpose!, BandIndex
+using Test: GenericArray
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
 
@@ -648,6 +651,16 @@ end
             @test_throws "cannot set index in the upper triangular part" copyto!(A, B)
         end
     end
+
+    @testset "partly initialized unit triangular" begin
+        for T in (UnitUpperTriangular, UnitLowerTriangular)
+            isupper = T == UnitUpperTriangular
+            M = Matrix{BigFloat}(undef, 2, 2)
+            M[1+!isupper,1+isupper] = 3
+            U = T(GenericArray(M))
+            @test copyto!(similar(M), U) == U
+        end
+    end
 end
 
 @testset "getindex with Integers" begin
@@ -923,4 +936,14 @@ end
     end
 end
 
+@testset "block unit triangular scaling" begin
+    m = SizedArrays.SizedArray{(2,2)}([1 2; 3 4])
+    U = UnitUpperTriangular(fill(m, 4, 4))
+    M = Matrix{eltype(U)}(U)
+    @test U/2 == M/2
+    @test 2\U == 2\M
+    @test U*2 == M*2
+    @test 2*U == 2*M
+end
+
 end # module TestTriangular
diff --git a/test/triangular2.jl b/test/triangular2.jl
index 4cad1373..f7e761af 100644
--- a/test/triangular2.jl
+++ b/test/triangular2.jl
@@ -2,6 +2,8 @@
 
 module TestTriangularReal
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasFloat, errorbounds, full!, transpose!
 
diff --git a/test/triangular3.jl b/test/triangular3.jl
index 82a2a147..a1990b57 100644
--- a/test/triangular3.jl
+++ b/test/triangular3.jl
@@ -2,6 +2,8 @@
 
 module TestTriangularComplex
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasFloat, errorbounds, full!, transpose!
 
diff --git a/test/tridiag.jl b/test/tridiag.jl
index 4b592a87..a7f0d7ad 100644
--- a/test/tridiag.jl
+++ b/test/tridiag.jl
@@ -2,6 +2,8 @@
 
 module TestTridiagonal
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
diff --git a/test/uniformscaling.jl b/test/uniformscaling.jl
index 10d427d1..21bd4b80 100644
--- a/test/uniformscaling.jl
+++ b/test/uniformscaling.jl
@@ -2,6 +2,8 @@
 
 module TestUniformscaling
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
diff --git a/test/unitful.jl b/test/unitful.jl
index cd4ade52..84e875e7 100644
--- a/test/unitful.jl
+++ b/test/unitful.jl
@@ -1,5 +1,7 @@
 module TestUnitfulLinAlg
 
+isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
+
 using Test, LinearAlgebra, Random
 
 Random.seed!(1234321)