From 9f63ee4f46b5ec9b2a680c055c912481b97e8d3f Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Sun, 1 Jun 2025 11:35:50 -0400
Subject: [PATCH 01/13] FillArrays extension

---
 Project.toml                         |  8 ++--
 ext/MatrixAlgebraKitFillArraysExt.jl | 71 ++++++++++++++++++++++++++++
 src/algorithms.jl                    |  1 -
 test/fillarrays.jl                   | 26 ++++++++++
 4 files changed, 102 insertions(+), 4 deletions(-)
 create mode 100644 ext/MatrixAlgebraKitFillArraysExt.jl
 create mode 100644 test/fillarrays.jl

diff --git a/Project.toml b/Project.toml
index 507d81df..28c8db66 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,21 +1,24 @@
 name = "MatrixAlgebraKit"
 uuid = "6c742aac-3347-4629-af66-fc926824e5e4"
 authors = ["Jutho <jutho.haegeman@ugent.be> and contributors"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 
 [extensions]
 MatrixAlgebraKitChainRulesCoreExt = "ChainRulesCore"
+MatrixAlgebraKitFillArraysExt = "FillArrays"
 
 [compat]
 Aqua = "0.6, 0.7, 0.8"
 ChainRulesCore = "1"
 ChainRulesTestUtils = "1"
+FillArrays = "1"
 JET = "0.9"
 LinearAlgebra = "1"
 SafeTestsets = "0.1"
@@ -36,5 +39,4 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore",
-        "ChainRulesTestUtils", "StableRNGs", "Zygote"]
+test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote"]
diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
new file mode 100644
index 00000000..c8603c00
--- /dev/null
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -0,0 +1,71 @@
+module MatrixAlgebraKitFillArraysExt
+
+using MatrixAlgebraKit
+using MatrixAlgebraKit: AbstractAlgorithm
+using FillArrays
+
+struct EyeAlgorithm <: AbstractAlgorithm end
+
+for f in [:eig_full,
+          :eigh_full,
+          :qr_compact,
+          :qr_full,
+          :left_polar,
+          :lq_compact,
+          :lq_full,
+          :right_polar,
+          :svd_compact,
+          :svd_full]
+    @eval begin
+        MatrixAlgebraKit.copy_input(::typeof($f), a::Eye) = a
+    end
+end
+
+for f in [:eig, :eigh, :lq, :qr, :polar, :svd]
+    ff = Symbol("default_", f, "_algorithm")
+    @eval begin
+        function MatrixAlgebraKit.$ff(a::Type{<:Eye}; kwargs...)
+            return EyeAlgorithm()
+        end
+    end
+end
+
+for f in [:eig_full!,
+          :eigh_full!,
+          :qr_compact!,
+          :qr_full!,
+          :left_polar!,
+          :lq_compact!,
+          :lq_full!,
+          :right_polar!]
+    @eval begin
+        nfactors(::typeof($f)) = 2
+    end
+end
+for f in [:svd_compact!, :svd_full!]
+    @eval begin
+        nfactors(::typeof($f)) = 3
+    end
+end
+
+for f in [:eig_full!,
+          :eigh_full!,
+          :qr_compact!,
+          :qr_full!,
+          :left_polar!,
+          :lq_compact!,
+          :lq_full!,
+          :right_polar!,
+          :svd_compact!,
+          :svd_full!]
+    @eval begin
+        function MatrixAlgebraKit.initialize_output(::typeof($f), a::Eye, alg::EyeAlgorithm)
+            return ntuple(_ -> a, nfactors($f))
+        end
+        function MatrixAlgebraKit.$f(a::Eye, F, alg::EyeAlgorithm; kwargs...)
+            return ntuple(_ -> a, nfactors($f))
+        end
+    end
+end
+
+end
diff --git a/src/algorithms.jl b/src/algorithms.jl
index f559b42e..6f9a4d49 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -99,7 +99,6 @@ function select_algorithm(f::F, ::Type{A}, alg::Alg=nothing; kwargs...) where {F
     throw(ArgumentError("Unknown alg $alg"))
 end
 
-
 @doc """
     MatrixAlgebraKit.default_algorithm(f, A; kwargs...)
     MatrixAlgebraKit.default_algorithm(f, ::Type{TA}; kwargs...) where {TA}
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
new file mode 100644
index 00000000..2b5c4b31
--- /dev/null
+++ b/test/fillarrays.jl
@@ -0,0 +1,26 @@
+using MatrixAlgebraKit
+using Test
+using TestExtras
+using FillArrays
+
+@testset "Eye" begin
+    for f in [:eig_full!,
+              :eigh_full!,
+              # TODO: Reenable once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32
+              # is merged.
+              # :qr_compact!,
+              # :qr_full!,
+              :left_polar!,
+              :lq_compact!,
+              :lq_full!,
+              :right_polar!,
+              :svd_compact!,
+              :svd_full!]
+        @eval begin
+            A = Eye(3)
+            F = @constinferred $f(A)
+            @test A === prod(F)
+            @test all(x -> x === A, F)
+        end
+    end
+end

From 28f627cb16692de60f5cb5b6f23c61a9709183be Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Sun, 1 Jun 2025 12:01:36 -0400
Subject: [PATCH 02/13] Support for Zeros

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 32 ++++++++++++++++++++++++++--
 test/fillarrays.jl                   | 13 +++++++----
 2 files changed, 39 insertions(+), 6 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index c8603c00..808e1496 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -1,10 +1,13 @@
 module MatrixAlgebraKitFillArraysExt
 
+using LinearAlgebra
 using MatrixAlgebraKit
 using MatrixAlgebraKit: AbstractAlgorithm
 using FillArrays
+using FillArrays: AbstractZerosMatrix
 
 struct EyeAlgorithm <: AbstractAlgorithm end
+struct ZerosAlgorithm <: AbstractAlgorithm end
 
 for f in [:eig_full,
           :eigh_full,
@@ -18,6 +21,8 @@ for f in [:eig_full,
           :svd_full]
     @eval begin
         MatrixAlgebraKit.copy_input(::typeof($f), a::Eye) = a
+
+        MatrixAlgebraKit.copy_input(::typeof($f), a::AbstractZerosMatrix) = a
     end
 end
 
@@ -27,6 +32,10 @@ for f in [:eig, :eigh, :lq, :qr, :polar, :svd]
         function MatrixAlgebraKit.$ff(a::Type{<:Eye}; kwargs...)
             return EyeAlgorithm()
         end
+
+        function MatrixAlgebraKit.$ff(a::Type{<:AbstractZerosMatrix}; kwargs...)
+            return ZerosAlgorithm()
+        end
     end
 end
 
@@ -59,12 +68,31 @@ for f in [:eig_full!,
           :svd_compact!,
           :svd_full!]
     @eval begin
-        function MatrixAlgebraKit.initialize_output(::typeof($f), a::Eye, alg::EyeAlgorithm)
-            return ntuple(_ -> a, nfactors($f))
+        function MatrixAlgebraKit.initialize_output(::typeof($f), a::Eye,
+                                                    alg::EyeAlgorithm)
+            return nothing
+        end
+        function MatrixAlgebraKit.check_input(::typeof($f), A::Eye, F)
+            LinearAlgebra.checksquare(A)
+            return nothing
         end
+
         function MatrixAlgebraKit.$f(a::Eye, F, alg::EyeAlgorithm; kwargs...)
             return ntuple(_ -> a, nfactors($f))
         end
+
+        function MatrixAlgebraKit.initialize_output(::typeof($f), a::AbstractZerosMatrix,
+                                                    alg::ZerosAlgorithm)
+            return nothing
+        end
+        function MatrixAlgebraKit.check_input(::typeof($f), A::AbstractZerosMatrix, F)
+            LinearAlgebra.checksquare(A)
+            return nothing
+        end
+        function MatrixAlgebraKit.$f(a::AbstractZerosMatrix, F, alg::ZerosAlgorithm;
+                                     kwargs...)
+            return ntuple(_ -> a, nfactors($f))
+        end
     end
 end
 
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index 2b5c4b31..ba24b728 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -6,10 +6,10 @@ using FillArrays
 @testset "Eye" begin
     for f in [:eig_full!,
               :eigh_full!,
-              # TODO: Reenable once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32
-              # is merged.
-              # :qr_compact!,
-              # :qr_full!,
+     # TODO: Reenable once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32
+     # is merged.
+     # :qr_compact!,
+     # :qr_full!,
               :left_polar!,
               :lq_compact!,
               :lq_full!,
@@ -21,6 +21,11 @@ using FillArrays
             F = @constinferred $f(A)
             @test A === prod(F)
             @test all(x -> x === A, F)
+
+            A = Zeros(3, 3)
+            F = @constinferred $f(A)
+            @test A === prod(F)
+            @test all(x -> x === A, F)
         end
     end
 end

From 1701e164df5704fa55732f326c342efc4e7a5e17 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Sun, 1 Jun 2025 19:49:55 -0400
Subject: [PATCH 03/13] More functionality and tests

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 208 ++++++++++++++++++++-------
 test/fillarrays.jl                   | 194 ++++++++++++++++++++++---
 2 files changed, 332 insertions(+), 70 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 808e1496..4bab590d 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -2,13 +2,21 @@ module MatrixAlgebraKitFillArraysExt
 
 using LinearAlgebra
 using MatrixAlgebraKit
-using MatrixAlgebraKit: AbstractAlgorithm
+using MatrixAlgebraKit: AbstractAlgorithm, check_input
 using FillArrays
 using FillArrays: AbstractZerosMatrix
 
-struct EyeAlgorithm <: AbstractAlgorithm end
 struct ZerosAlgorithm <: AbstractAlgorithm end
 
+for f in [:eig, :eigh, :lq, :qr, :svd]
+    ff = Symbol("default_", f, "_algorithm")
+    @eval begin
+        function MatrixAlgebraKit.$ff(::Type{<:AbstractZerosMatrix}; kwargs...)
+            return ZerosAlgorithm()
+        end
+    end
+end
+
 for f in [:eig_full,
           :eigh_full,
           :qr_compact,
@@ -19,81 +27,181 @@ for f in [:eig_full,
           :right_polar,
           :svd_compact,
           :svd_full]
+    f! = Symbol(f, "!")
     @eval begin
-        MatrixAlgebraKit.copy_input(::typeof($f), a::Eye) = a
-
-        MatrixAlgebraKit.copy_input(::typeof($f), a::AbstractZerosMatrix) = a
+        MatrixAlgebraKit.copy_input(::typeof($f), A::AbstractZerosMatrix) = A
+        function MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractZerosMatrix,
+                                                    alg::ZerosAlgorithm)
+            return nothing
+        end
     end
 end
 
-for f in [:eig, :eigh, :lq, :qr, :polar, :svd]
-    ff = Symbol("default_", f, "_algorithm")
+for f in [:eig_full!,
+          :eigh_full!]
     @eval begin
-        function MatrixAlgebraKit.$ff(a::Type{<:Eye}; kwargs...)
-            return EyeAlgorithm()
+        function MatrixAlgebraKit.check_input(::typeof($f), A::AbstractZerosMatrix, F)
+            LinearAlgebra.checksquare(A)
+            return nothing
         end
-
-        function MatrixAlgebraKit.$ff(a::Type{<:AbstractZerosMatrix}; kwargs...)
-            return ZerosAlgorithm()
+        function MatrixAlgebraKit.$f(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm;
+                                     kwargs...)
+            check_input($f, A, F)
+            return (A, Eye(axes(A)))
         end
     end
 end
 
-for f in [:eig_full!,
-          :eigh_full!,
-          :qr_compact!,
-          :qr_full!,
-          :left_polar!,
-          :lq_compact!,
-          :lq_full!,
-          :right_polar!]
-    @eval begin
-        nfactors(::typeof($f)) = 2
+function MatrixAlgebraKit.qr_compact!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    m, n = size(A)
+    ax = axes(A)
+    if m > n
+        r_ax = (ax[2], ax[2])
+        return (Eye(ax), Zeros(r_ax))
+    else
+        q_ax = (ax[1], ax[1])
+        return (Eye(q_ax), Zeros(ax))
     end
 end
-for f in [:svd_compact!, :svd_full!]
+
+function MatrixAlgebraKit.qr_full!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    ax = axes(A)
+    q_ax = (ax[1], ax[1])
+    return (Eye(q_ax), Zeros(ax))
+end
+
+function MatrixAlgebraKit.lq_compact!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    m, n = size(A)
+    ax = axes(A)
+    if m < n
+        l_ax = (ax[1], ax[1])
+        return (Zeros(l_ax), Eye(ax))
+    else
+        q_ax = (ax[2], ax[2])
+        return (Zeros(ax), Eye(q_ax))
+    end
+end
+
+function MatrixAlgebraKit.lq_full!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    ax = axes(A)
+    q_ax = (ax[2], ax[2])
+    return (Zeros(ax), Eye(q_ax))
+end
+
+function MatrixAlgebraKit.svd_compact!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    m, n = size(A)
+    ax = axes(A)
+    if m > n
+        s_ax = (ax[2], ax[2])
+        return (Eye(ax), Zeros(s_ax), Eye(s_ax))
+    else
+        s_ax = (ax[1], ax[1])
+        return (Eye(s_ax), Zeros(s_ax), Eye(ax))
+    end
+end
+
+function MatrixAlgebraKit.svd_full!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    ax = axes(A)
+    return (Eye((ax[1], ax[1])), Zeros(ax), Eye((ax[2], ax[2])))
+end
+
+struct EyeAlgorithm <: AbstractAlgorithm end
+
+for f in [:eig, :eigh, :lq, :qr, :polar, :svd]
+    ff = Symbol("default_", f, "_algorithm")
     @eval begin
-        nfactors(::typeof($f)) = 3
+        function MatrixAlgebraKit.$ff(A::Type{<:Eye}; kwargs...)
+            return EyeAlgorithm()
+        end
     end
 end
 
-for f in [:eig_full!,
-          :eigh_full!,
-          :qr_compact!,
-          :qr_full!,
-          :left_polar!,
-          :lq_compact!,
-          :lq_full!,
-          :right_polar!,
-          :svd_compact!,
-          :svd_full!]
+for f in [:eig_full,
+          :eigh_full,
+          :qr_compact,
+          :qr_full,
+          :lq_compact,
+          :lq_full,
+          :left_polar,
+          :right_polar,
+          :svd_compact,
+          :svd_full]
+    f! = Symbol(f, "!")
     @eval begin
-        function MatrixAlgebraKit.initialize_output(::typeof($f), a::Eye,
+        MatrixAlgebraKit.copy_input(::typeof($f), A::Eye) = A
+        function MatrixAlgebraKit.initialize_output(::typeof($f!), A::Eye,
                                                     alg::EyeAlgorithm)
             return nothing
         end
-        function MatrixAlgebraKit.check_input(::typeof($f), A::Eye, F)
-            LinearAlgebra.checksquare(A)
-            return nothing
-        end
-
-        function MatrixAlgebraKit.$f(a::Eye, F, alg::EyeAlgorithm; kwargs...)
-            return ntuple(_ -> a, nfactors($f))
-        end
+    end
+end
 
-        function MatrixAlgebraKit.initialize_output(::typeof($f), a::AbstractZerosMatrix,
-                                                    alg::ZerosAlgorithm)
-            return nothing
-        end
-        function MatrixAlgebraKit.check_input(::typeof($f), A::AbstractZerosMatrix, F)
+for f in [:eig_full!,
+          :eigh_full!]
+    @eval begin
+        function MatrixAlgebraKit.check_input(::typeof($f), A::Eye, F)
             LinearAlgebra.checksquare(A)
             return nothing
         end
-        function MatrixAlgebraKit.$f(a::AbstractZerosMatrix, F, alg::ZerosAlgorithm;
+        function MatrixAlgebraKit.$f(A::Eye, F, alg::EyeAlgorithm;
                                      kwargs...)
-            return ntuple(_ -> a, nfactors($f))
+            check_input($f, A, F)
+            return (A, A)
         end
     end
 end
 
+function MatrixAlgebraKit.qr_compact!(A::Eye, F, alg::EyeAlgorithm)
+    m, n = size(A)
+    ax = axes(A)
+    if m > n
+        r_ax = (ax[2], ax[2])
+        return (Eye(ax), Eye(r_ax))
+    else
+        q_ax = (ax[1], ax[1])
+        return (Eye(q_ax), Eye(ax))
+    end
+end
+
+function MatrixAlgebraKit.qr_full!(A::Eye, F, alg::EyeAlgorithm)
+    ax = axes(A)
+    q_ax = (ax[1], ax[1])
+    return (Eye(q_ax), A)
+end
+
+function MatrixAlgebraKit.lq_compact!(A::Eye, F, alg::EyeAlgorithm)
+    m, n = size(A)
+    ax = axes(A)
+    if m < n
+        l_ax = (ax[1], ax[1])
+        return (Eye(l_ax), Eye(ax))
+    else
+        q_ax = (ax[2], ax[2])
+        return (Eye(ax), Eye(q_ax))
+    end
+end
+
+function MatrixAlgebraKit.lq_full!(A::Eye, F, alg::EyeAlgorithm)
+    ax = axes(A)
+    q_ax = (ax[2], ax[2])
+    return (A, Eye(q_ax))
+end
+
+function MatrixAlgebraKit.svd_compact!(A::Eye, F, alg::EyeAlgorithm)
+    m, n = size(A)
+    ax = axes(A)
+    if m > n
+        s_ax = (ax[2], ax[2])
+        return (Eye(ax), Eye(s_ax), Eye(s_ax))
+    else
+        s_ax = (ax[1], ax[1])
+        return (Eye(s_ax), Eye(s_ax), Eye(ax))
+    end
+end
+
+function MatrixAlgebraKit.svd_full!(A::Eye, F, alg::EyeAlgorithm)
+    ax = axes(A)
+    return (Eye((ax[1], ax[1])), A, Eye((ax[2], ax[2])))
+end
+
 end
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index ba24b728..3a75262b 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -1,31 +1,185 @@
 using MatrixAlgebraKit
+using LinearAlgebra
 using Test
 using TestExtras
 using FillArrays
 
+@testset "Zeros" begin
+    for f in [:eig_full, :eigh_full]
+        @eval begin
+            A = Zeros(3, 3)
+            D, V = @constinferred $f(A)
+            @test A * V == D * V
+            @test size(D) == size(A)
+            @test size(V) == size(A)
+            @test iszero(D)
+            @test D isa Zeros
+            @test V == I
+            @test V isa Eye
+        end
+    end
+
+    # for f in [:qr_compact, :left_polar]
+    #     @eval begin
+    #         A = Zeros(4, 3)
+    #         Q, R = @constinferred $f(A)
+    #         @test A == Q * R
+    #         @test size(Q) == (4, 3)
+    #         @test size(R) == (3, 3)
+    #         @test Q == Matrix(I, (4, 3))
+    #         @test Q isa Eye
+    #         @test iszero(R)
+    #         @test R isa Zeros
+    #     end
+    # end
+
+    # A = Zeros(4, 3)
+    # Q, R = @constinferred qr_full(A)
+    # @test A == Q * R
+    # @test size(Q) == (4, 4)
+    # @test size(R) == (4, 3)
+    # @test Q == I
+    # @test Q isa Eye
+    # @test iszero(R)
+    # @test R isa Zeros
+
+    for f in [:lq_compact] # :right_polar]
+        @eval begin
+            A = Zeros(3, 4)
+            L, Q = @constinferred $f(A)
+            @test A == L * Q
+            @test size(L) == (3, 3)
+            @test size(Q) == (3, 4)
+            @test iszero(L)
+            @test L isa Zeros
+            @test Q == Matrix(I, (3, 4))
+            @test Q isa Eye
+        end
+    end
+
+    A = Zeros(3, 4)
+    L, Q = @constinferred lq_full(A)
+    @test A == L * Q
+    @test size(L) == (3, 4)
+    @test size(Q) == (4, 4)
+    @test iszero(L)
+    @test L isa Zeros
+    @test Q == I
+    @test Q isa Eye
+
+    A = Zeros(3, 4)
+    U, S, V = @constinferred svd_compact(A)
+    @test U * S * V == A
+    @test size(U) == (3, 3)
+    @test size(S) == (3, 3)
+    @test size(V) == (3, 4)
+    @test iszero(S)
+    @test S isa Zeros
+    @test U == I
+    @test U isa Eye
+    @test V == Matrix(I, (3, 4))
+    @test V isa Eye
+
+    A = Zeros(3, 4)
+    U, S, V = @constinferred svd_full(A)
+    @test U * S * V == A
+    @test size(U) == (3, 3)
+    @test size(S) == (3, 4)
+    @test size(V) == (4, 4)
+    @test iszero(S)
+    @test S isa Zeros
+    @test U == I
+    @test U isa Eye
+    @test V == I
+    @test V isa Eye
+end
+
 @testset "Eye" begin
-    for f in [:eig_full!,
-              :eigh_full!,
-     # TODO: Reenable once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32
-     # is merged.
-     # :qr_compact!,
-     # :qr_full!,
-              :left_polar!,
-              :lq_compact!,
-              :lq_full!,
-              :right_polar!,
-              :svd_compact!,
-              :svd_full!]
+    for f in [:eig_full, :eigh_full]
         @eval begin
-            A = Eye(3)
-            F = @constinferred $f(A)
-            @test A === prod(F)
-            @test all(x -> x === A, F)
+            A = Eye(3, 3)
+            D, V = @constinferred $f(A)
+            @test A * V == D * V
+            @test size(D) == size(A)
+            @test size(V) == size(A)
+            @test V == I
+            @test D isa Eye
+            @test V == I
+            @test V isa Eye
+        end
+    end
 
-            A = Zeros(3, 3)
-            F = @constinferred $f(A)
-            @test A === prod(F)
-            @test all(x -> x === A, F)
+    # for f in [:qr_compact, :left_polar]
+    #     @eval begin
+    #         A = Eye(4, 3)
+    #         Q, R = @constinferred $f(A)
+    #         @test A == Q * R
+    #         @test size(Q) == (4, 3)
+    #         @test size(R) == (3, 3)
+    #         @test Q == Matrix(I, (4, 3))
+    #         @test Q isa Eye
+    #         @test R == I
+    #         @test R isa Eye
+    #     end
+    # end
+
+    # A = Eye(4, 3)
+    # Q, R = @constinferred qr_full(A)
+    # @test A == Q * R
+    # @test size(Q) == (4, 4)
+    # @test size(R) == (4, 3)
+    # @test Q == I
+    # @test Q isa Eye
+    # @test R == I
+    # @test R isa Eye
+
+    for f in [:lq_compact] # :right_polar]
+        @eval begin
+            A = Eye(3, 4)
+            L, Q = @constinferred $f(A)
+            @test A == L * Q
+            @test size(L) == (3, 3)
+            @test size(Q) == (3, 4)
+            @test L == I
+            @test L isa Eye
+            @test Q == Matrix(I, (3, 4))
+            @test Q isa Eye
         end
     end
+
+    A = Eye(3, 4)
+    L, Q = @constinferred lq_full(A)
+    @test A == L * Q
+    @test size(L) == (3, 4)
+    @test size(Q) == (4, 4)
+    @test L == Matrix(I, (3, 4))
+    @test L isa Eye
+    @test Q == I
+    @test Q isa Eye
+
+    A = Eye(3, 4)
+    U, S, V = @constinferred svd_compact(A)
+    @test U * S * V == A
+    @test size(U) == (3, 3)
+    @test size(S) == (3, 3)
+    @test size(V) == (3, 4)
+    @test S == I
+    @test S isa Eye
+    @test U == I
+    @test U isa Eye
+    @test V == Matrix(I, (3, 4))
+    @test V isa Eye
+
+    A = Eye(3, 4)
+    U, S, V = @constinferred svd_full(A)
+    @test U * S * V == A
+    @test size(U) == (3, 3)
+    @test size(S) == (3, 4)
+    @test size(V) == (4, 4)
+    @test S == Matrix(I, (3, 4))
+    @test S isa Eye
+    @test U == I
+    @test U isa Eye
+    @test V == I
+    @test V isa Eye
 end

From 78b841cee49d0187254e6da91b53aa2bc7d8e592 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Sun, 1 Jun 2025 20:05:00 -0400
Subject: [PATCH 04/13] Start adding support for eig_vals

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 46 +++++++++++++++++++++++-----
 test/fillarrays.jl                   | 20 ++++++++++++
 2 files changed, 59 insertions(+), 7 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 4bab590d..67566688 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -2,7 +2,7 @@ module MatrixAlgebraKitFillArraysExt
 
 using LinearAlgebra
 using MatrixAlgebraKit
-using MatrixAlgebraKit: AbstractAlgorithm, check_input
+using MatrixAlgebraKit: AbstractAlgorithm, check_input, diagview
 using FillArrays
 using FillArrays: AbstractZerosMatrix
 
@@ -19,6 +19,8 @@ end
 
 for f in [:eig_full,
           :eigh_full,
+          :eig_vals,
+          :eigh_vals,
           :qr_compact,
           :qr_full,
           :left_polar,
@@ -26,7 +28,8 @@ for f in [:eig_full,
           :lq_full,
           :right_polar,
           :svd_compact,
-          :svd_full]
+          :svd_full,
+          :svd_vals]
     f! = Symbol(f, "!")
     @eval begin
         MatrixAlgebraKit.copy_input(::typeof($f), A::AbstractZerosMatrix) = A
@@ -37,8 +40,7 @@ for f in [:eig_full,
     end
 end
 
-for f in [:eig_full!,
-          :eigh_full!]
+for f in [:eig_full!, :eigh_full!]
     @eval begin
         function MatrixAlgebraKit.check_input(::typeof($f), A::AbstractZerosMatrix, F)
             LinearAlgebra.checksquare(A)
@@ -52,6 +54,20 @@ for f in [:eig_full!,
     end
 end
 
+for f in [:eig_vals!, :eigh_vals!]
+    @eval begin
+        function MatrixAlgebraKit.check_input(::typeof($f), A::AbstractZerosMatrix, F)
+            LinearAlgebra.checksquare(A)
+            return nothing
+        end
+        function MatrixAlgebraKit.$f(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm;
+                                     kwargs...)
+            check_input($f, A, F)
+            return Zeros(axes(A, 1))
+        end
+    end
+end
+
 function MatrixAlgebraKit.qr_compact!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
     m, n = size(A)
     ax = axes(A)
@@ -118,6 +134,8 @@ end
 
 for f in [:eig_full,
           :eigh_full,
+          :eig_vals,
+          :eigh_vals,
           :qr_compact,
           :qr_full,
           :lq_compact,
@@ -125,7 +143,8 @@ for f in [:eig_full,
           :left_polar,
           :right_polar,
           :svd_compact,
-          :svd_full]
+          :svd_full,
+          :svd_vals]
     f! = Symbol(f, "!")
     @eval begin
         MatrixAlgebraKit.copy_input(::typeof($f), A::Eye) = A
@@ -136,8 +155,7 @@ for f in [:eig_full,
     end
 end
 
-for f in [:eig_full!,
-          :eigh_full!]
+for f in [:eig_full!, :eigh_full!]
     @eval begin
         function MatrixAlgebraKit.check_input(::typeof($f), A::Eye, F)
             LinearAlgebra.checksquare(A)
@@ -151,6 +169,20 @@ for f in [:eig_full!,
     end
 end
 
+for f in [:eig_vals!, :eigh_vals!]
+    @eval begin
+        function MatrixAlgebraKit.check_input(::typeof($f), A::Eye, F)
+            LinearAlgebra.checksquare(A)
+            return nothing
+        end
+        function MatrixAlgebraKit.$f(A::Eye, F, alg::EyeAlgorithm;
+                                     kwargs...)
+            check_input($f, A, F)
+            return diagview(A)
+        end
+    end
+end
+
 function MatrixAlgebraKit.qr_compact!(A::Eye, F, alg::EyeAlgorithm)
     m, n = size(A)
     ax = axes(A)
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index 3a75262b..d08132d0 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -19,6 +19,16 @@ using FillArrays
         end
     end
 
+    for f in [:eig_vals, :eigh_vals]
+        @eval begin
+            A = Zeros(3, 3)
+            D = @constinferred $f(A)
+            @test size(D) == (size(A, 1),)
+            @test iszero(D)
+            @test D isa Zeros
+        end
+    end
+
     # for f in [:qr_compact, :left_polar]
     #     @eval begin
     #         A = Zeros(4, 3)
@@ -109,6 +119,16 @@ end
         end
     end
 
+    for f in [:eig_vals, :eigh_vals]
+        @eval begin
+            A = Eye(3, 3)
+            D = @constinferred $f(A)
+            @test size(D) == (size(A, 1),)
+            @test all(isone, D)
+            @test D isa Ones
+        end
+    end
+
     # for f in [:qr_compact, :left_polar]
     #     @eval begin
     #         A = Eye(4, 3)

From c963e5fb835a4b931b9a805c489cc7bc0813e459 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Sun, 1 Jun 2025 20:46:59 -0400
Subject: [PATCH 05/13] svd_vals

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 16 ++++++++++++++--
 test/fillarrays.jl                   | 12 ++++++++++++
 2 files changed, 26 insertions(+), 2 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 67566688..d2165aba 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -4,7 +4,11 @@ using LinearAlgebra
 using MatrixAlgebraKit
 using MatrixAlgebraKit: AbstractAlgorithm, check_input, diagview
 using FillArrays
-using FillArrays: AbstractZerosMatrix
+using FillArrays: AbstractZerosMatrix, OnesVector, RectDiagonal
+
+function MatrixAlgebraKit.diagview(A::RectDiagonal{<:Any,<:OnesVector})
+    return A.diag
+end
 
 struct ZerosAlgorithm <: AbstractAlgorithm end
 
@@ -63,7 +67,7 @@ for f in [:eig_vals!, :eigh_vals!]
         function MatrixAlgebraKit.$f(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm;
                                      kwargs...)
             check_input($f, A, F)
-            return Zeros(axes(A, 1))
+            return diagview(A)
         end
     end
 end
@@ -121,6 +125,10 @@ function MatrixAlgebraKit.svd_full!(A::AbstractZerosMatrix, F, alg::ZerosAlgorit
     return (Eye((ax[1], ax[1])), Zeros(ax), Eye((ax[2], ax[2])))
 end
 
+function MatrixAlgebraKit.svd_vals!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    return diagview(A)
+end
+
 struct EyeAlgorithm <: AbstractAlgorithm end
 
 for f in [:eig, :eigh, :lq, :qr, :polar, :svd]
@@ -236,4 +244,8 @@ function MatrixAlgebraKit.svd_full!(A::Eye, F, alg::EyeAlgorithm)
     return (Eye((ax[1], ax[1])), A, Eye((ax[2], ax[2])))
 end
 
+function MatrixAlgebraKit.svd_vals!(A::Eye, F, alg::EyeAlgorithm)
+    return diagview(A)
+end
+
 end
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index d08132d0..a5d1e757 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -102,6 +102,12 @@ using FillArrays
     @test U isa Eye
     @test V == I
     @test V isa Eye
+
+    A = Zeros(3, 4)
+    D = @constinferred svd_vals(A)
+    @test size(D) == (minimum(size(A)),)
+    @test iszero(D)
+    @test D isa Zeros
 end
 
 @testset "Eye" begin
@@ -202,4 +208,10 @@ end
     @test U isa Eye
     @test V == I
     @test V isa Eye
+
+    A = Eye(3, 4)
+    D = @constinferred svd_vals(A)
+    @test size(D) == (minimum(size(A)),)
+    @test all(isone, D)
+    @test D isa Ones
 end

From aa523c96e79dc04ff1b61a94e2c8a2ac9e52f275 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Sun, 1 Jun 2025 21:18:30 -0400
Subject: [PATCH 06/13] Specialize on square Eye

---
 ext/MatrixAlgebraKitFillArraysExt.jl |  20 ++-
 test/fillarrays.jl                   | 223 ++++++++++++++++++---------
 2 files changed, 172 insertions(+), 71 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index d2165aba..9cc7dee8 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -4,7 +4,7 @@ using LinearAlgebra
 using MatrixAlgebraKit
 using MatrixAlgebraKit: AbstractAlgorithm, check_input, diagview
 using FillArrays
-using FillArrays: AbstractZerosMatrix, OnesVector, RectDiagonal
+using FillArrays: AbstractZerosMatrix, OnesVector, RectDiagonal, SquareEye
 
 function MatrixAlgebraKit.diagview(A::RectDiagonal{<:Any,<:OnesVector})
     return A.diag
@@ -202,12 +202,18 @@ function MatrixAlgebraKit.qr_compact!(A::Eye, F, alg::EyeAlgorithm)
         return (Eye(q_ax), Eye(ax))
     end
 end
+function MatrixAlgebraKit.qr_compact!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A)
+end
 
 function MatrixAlgebraKit.qr_full!(A::Eye, F, alg::EyeAlgorithm)
     ax = axes(A)
     q_ax = (ax[1], ax[1])
     return (Eye(q_ax), A)
 end
+function MatrixAlgebraKit.qr_full!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A)
+end
 
 function MatrixAlgebraKit.lq_compact!(A::Eye, F, alg::EyeAlgorithm)
     m, n = size(A)
@@ -220,12 +226,18 @@ function MatrixAlgebraKit.lq_compact!(A::Eye, F, alg::EyeAlgorithm)
         return (Eye(ax), Eye(q_ax))
     end
 end
+function MatrixAlgebraKit.lq_compact!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A)
+end
 
 function MatrixAlgebraKit.lq_full!(A::Eye, F, alg::EyeAlgorithm)
     ax = axes(A)
     q_ax = (ax[2], ax[2])
     return (A, Eye(q_ax))
 end
+function MatrixAlgebraKit.lq_full!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A)
+end
 
 function MatrixAlgebraKit.svd_compact!(A::Eye, F, alg::EyeAlgorithm)
     m, n = size(A)
@@ -238,11 +250,17 @@ function MatrixAlgebraKit.svd_compact!(A::Eye, F, alg::EyeAlgorithm)
         return (Eye(s_ax), Eye(s_ax), Eye(ax))
     end
 end
+function MatrixAlgebraKit.svd_compact!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A, A)
+end
 
 function MatrixAlgebraKit.svd_full!(A::Eye, F, alg::EyeAlgorithm)
     ax = axes(A)
     return (Eye((ax[1], ax[1])), A, Eye((ax[2], ax[2])))
 end
+function MatrixAlgebraKit.svd_full!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A, A)
+end
 
 function MatrixAlgebraKit.svd_vals!(A::Eye, F, alg::EyeAlgorithm)
     return diagview(A)
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index a5d1e757..f1d8dc14 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -29,23 +29,25 @@ using FillArrays
         end
     end
 
-    # for f in [:qr_compact, :left_polar]
-    #     @eval begin
-    #         A = Zeros(4, 3)
-    #         Q, R = @constinferred $f(A)
-    #         @test A == Q * R
-    #         @test size(Q) == (4, 3)
-    #         @test size(R) == (3, 3)
-    #         @test Q == Matrix(I, (4, 3))
-    #         @test Q isa Eye
-    #         @test iszero(R)
-    #         @test R isa Zeros
-    #     end
-    # end
+    # A = Zeros(4, 3)
+    # Q, R = @constinferred qr_compact(A)
+    # @test Q * R == A
+    # @test size(Q) == (4, 3)
+    # @test size(R) == (3, 3)
+    # @test Q == Matrix(I, (4, 3))
+    # @test Q isa Eye
+    # @test iszero(R)
+    # @test R isa Zeros
+
+    # A = Zeros(3)
+    # Q, R = @constinferred qr_compact(A)
+    # @test Q * R == A
+    # @test Q === A
+    # @test R === A
 
     # A = Zeros(4, 3)
     # Q, R = @constinferred qr_full(A)
-    # @test A == Q * R
+    # @test Q * R == A
     # @test size(Q) == (4, 4)
     # @test size(R) == (4, 3)
     # @test Q == I
@@ -53,23 +55,49 @@ using FillArrays
     # @test iszero(R)
     # @test R isa Zeros
 
-    for f in [:lq_compact] # :right_polar]
-        @eval begin
-            A = Zeros(3, 4)
-            L, Q = @constinferred $f(A)
-            @test A == L * Q
-            @test size(L) == (3, 3)
-            @test size(Q) == (3, 4)
-            @test iszero(L)
-            @test L isa Zeros
-            @test Q == Matrix(I, (3, 4))
-            @test Q isa Eye
-        end
-    end
+    # A = Zeros(3)
+    # Q, R = @constinferred qr_full(A)
+    # @test Q * R == A
+    # @test Q === A
+    # @test R === A
+
+    # A = Zeros(4, 3)
+    # Q, R = @constinferred left_polar(A)
+    # @test Q * R == A
+    # @test size(Q) == (4, 3)
+    # @test size(R) == (3, 3)
+    # @test Q == Matrix(I, (4, 3))
+    # @test Q isa Eye
+    # @test iszero(R)
+    # @test R isa Zeros
+
+    # A = Zeros(3)
+    # Q, R = @constinferred left_polar(A)
+    # @test Q * R == A
+    # @test Q === A
+    # @test R === A
+
+    A = Zeros(3, 4)
+    L, Q = @constinferred lq_compact(A)
+    @test L * Q == A
+    @test size(L) == (3, 3)
+    @test size(Q) == (3, 4)
+    @test iszero(L)
+    @test L isa Zeros
+    @test Q == Matrix(I, (3, 4))
+    @test Q isa Eye
+
+    A = Zeros(3, 4)
+    L, Q = @constinferred lq_compact(A)
+    @test L * Q == A
+    @test L == Zeros(3, 3)
+    @test L isa Zeros
+    @test Q == Eye(3, 4)
+    @test Q isa Eye
 
     A = Zeros(3, 4)
     L, Q = @constinferred lq_full(A)
-    @test A == L * Q
+    @test L * Q == A
     @test size(L) == (3, 4)
     @test size(Q) == (4, 4)
     @test iszero(L)
@@ -77,6 +105,29 @@ using FillArrays
     @test Q == I
     @test Q isa Eye
 
+    A = Zeros(3, 4)
+    L, Q = @constinferred lq_full(A)
+    @test L * Q == A
+    @test L === A
+    @test Q == Eye(4)
+    @test Q isa Eye
+
+    # A = Zeros(3, 4)
+    # L, Q = @constinferred right_polar(A)
+    # @test L * Q == A
+    # @test size(L) == (3, 3)
+    # @test size(Q) == (3, 4)
+    # @test iszero(L)
+    # @test L isa Zeros
+    # @test Q == Matrix(I, (3, 4))
+    # @test Q isa Eye
+
+    # A = Zeros(3, 4)
+    # L, Q = @constinferred right_polar(A)
+    # @test L * Q == A
+    # @test L === A
+    # @test Q === A
+
     A = Zeros(3, 4)
     U, S, V = @constinferred svd_compact(A)
     @test U * S * V == A
@@ -113,45 +164,49 @@ end
 @testset "Eye" begin
     for f in [:eig_full, :eigh_full]
         @eval begin
-            A = Eye(3, 3)
-            D, V = @constinferred $f(A)
-            @test A * V == D * V
-            @test size(D) == size(A)
-            @test size(V) == size(A)
-            @test V == I
-            @test D isa Eye
-            @test V == I
-            @test V isa Eye
+            for A in (Eye(3), Eye(3, 3))
+                local D, V = @constinferred $f(A)
+                @test A * V == D * V
+                @test size(D) == size(A)
+                @test size(V) == size(A)
+                @test V == I
+                @test typeof(D) === typeof(A)
+                @test V == I
+                @test typeof(V) === typeof(A)
+            end
         end
     end
 
     for f in [:eig_vals, :eigh_vals]
         @eval begin
-            A = Eye(3, 3)
-            D = @constinferred $f(A)
-            @test size(D) == (size(A, 1),)
-            @test all(isone, D)
-            @test D isa Ones
+            for A in (Eye(3), Eye(3, 3))
+                local D = @constinferred $f(A)
+                @test size(D) == (size(A, 1),)
+                @test all(isone, D)
+                @test D isa Ones
+            end
         end
     end
 
-    # for f in [:qr_compact, :left_polar]
-    #     @eval begin
-    #         A = Eye(4, 3)
-    #         Q, R = @constinferred $f(A)
-    #         @test A == Q * R
-    #         @test size(Q) == (4, 3)
-    #         @test size(R) == (3, 3)
-    #         @test Q == Matrix(I, (4, 3))
-    #         @test Q isa Eye
-    #         @test R == I
-    #         @test R isa Eye
-    #     end
-    # end
+    # A = Eye(4, 3)
+    # Q, R = @constinferred qr_compact(A)
+    # @test Q * R == A
+    # @test size(Q) == (4, 3)
+    # @test size(R) == (3, 3)
+    # @test Q == Matrix(I, (4, 3))
+    # @test Q isa Eye
+    # @test R == I
+    # @test R isa Eye
+
+    # A = Eye(3)
+    # Q, R = @constinferred qr_compact(A)
+    # @test Q * R == A
+    # @test Q === A
+    # @test R === A
 
     # A = Eye(4, 3)
     # Q, R = @constinferred qr_full(A)
-    # @test A == Q * R
+    # @test Q * R == A
     # @test size(Q) == (4, 4)
     # @test size(R) == (4, 3)
     # @test Q == I
@@ -159,23 +214,31 @@ end
     # @test R == I
     # @test R isa Eye
 
-    for f in [:lq_compact] # :right_polar]
-        @eval begin
-            A = Eye(3, 4)
-            L, Q = @constinferred $f(A)
-            @test A == L * Q
-            @test size(L) == (3, 3)
-            @test size(Q) == (3, 4)
-            @test L == I
-            @test L isa Eye
-            @test Q == Matrix(I, (3, 4))
-            @test Q isa Eye
-        end
-    end
+    # A = Eye(3)
+    # Q, R = @constinferred qr_full(A)
+    # @test Q * R == A
+    # @test Q === A
+    # @test R === A
+
+    A = Eye(3, 4)
+    L, Q = @constinferred lq_compact(A)
+    @test L * Q == A
+    @test size(L) == (3, 3)
+    @test size(Q) == (3, 4)
+    @test L == I
+    @test L isa Eye
+    @test Q == Matrix(I, (3, 4))
+    @test Q isa Eye
+
+    A = Eye(3)
+    L, Q = @constinferred lq_compact(A)
+    @test L * Q == A
+    @test L === A
+    @test Q === A
 
     A = Eye(3, 4)
     L, Q = @constinferred lq_full(A)
-    @test A == L * Q
+    @test L * Q == A
     @test size(L) == (3, 4)
     @test size(Q) == (4, 4)
     @test L == Matrix(I, (3, 4))
@@ -183,6 +246,12 @@ end
     @test Q == I
     @test Q isa Eye
 
+    A = Eye(3)
+    L, Q = @constinferred lq_full(A)
+    @test L * Q == A
+    @test L === A
+    @test Q === A
+
     A = Eye(3, 4)
     U, S, V = @constinferred svd_compact(A)
     @test U * S * V == A
@@ -196,6 +265,13 @@ end
     @test V == Matrix(I, (3, 4))
     @test V isa Eye
 
+    A = Eye(3)
+    U, S, V = @constinferred svd_compact(A)
+    @test U * S * V == A
+    @test U === A
+    @test S === A
+    @test V === A
+
     A = Eye(3, 4)
     U, S, V = @constinferred svd_full(A)
     @test U * S * V == A
@@ -209,6 +285,13 @@ end
     @test V == I
     @test V isa Eye
 
+    A = Eye(3)
+    U, S, V = @constinferred svd_full(A)
+    @test U * S * V == A
+    @test U === A
+    @test S === A
+    @test V === A
+
     A = Eye(3, 4)
     D = @constinferred svd_vals(A)
     @test size(D) == (minimum(size(A)),)

From a2658b2c2dc554bfc197d493b15d87d2df5e2be5 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 2 Jun 2025 07:48:36 -0400
Subject: [PATCH 07/13] More tests

---
 ext/MatrixAlgebraKitFillArraysExt.jl |  2 +-
 test/fillarrays.jl                   | 56 ++++++++++++++++------------
 2 files changed, 33 insertions(+), 25 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 9cc7dee8..62f729bf 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -256,7 +256,7 @@ end
 
 function MatrixAlgebraKit.svd_full!(A::Eye, F, alg::EyeAlgorithm)
     ax = axes(A)
-    return (Eye((ax[1], ax[1])), A, Eye((ax[2], ax[2])))
+    return (Eye((ax[1],)), A, Eye((ax[2],)))
 end
 function MatrixAlgebraKit.svd_full!(A::SquareEye, F, alg::EyeAlgorithm)
     return (A, A, A)
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index f1d8dc14..f772d641 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -39,12 +39,6 @@ using FillArrays
     # @test iszero(R)
     # @test R isa Zeros
 
-    # A = Zeros(3)
-    # Q, R = @constinferred qr_compact(A)
-    # @test Q * R == A
-    # @test Q === A
-    # @test R === A
-
     # A = Zeros(4, 3)
     # Q, R = @constinferred qr_full(A)
     # @test Q * R == A
@@ -55,12 +49,6 @@ using FillArrays
     # @test iszero(R)
     # @test R isa Zeros
 
-    # A = Zeros(3)
-    # Q, R = @constinferred qr_full(A)
-    # @test Q * R == A
-    # @test Q === A
-    # @test R === A
-
     # A = Zeros(4, 3)
     # Q, R = @constinferred left_polar(A)
     # @test Q * R == A
@@ -71,12 +59,6 @@ using FillArrays
     # @test iszero(R)
     # @test R isa Zeros
 
-    # A = Zeros(3)
-    # Q, R = @constinferred left_polar(A)
-    # @test Q * R == A
-    # @test Q === A
-    # @test R === A
-
     A = Zeros(3, 4)
     L, Q = @constinferred lq_compact(A)
     @test L * Q == A
@@ -122,12 +104,6 @@ using FillArrays
     # @test Q == Matrix(I, (3, 4))
     # @test Q isa Eye
 
-    # A = Zeros(3, 4)
-    # L, Q = @constinferred right_polar(A)
-    # @test L * Q == A
-    # @test L === A
-    # @test Q === A
-
     A = Zeros(3, 4)
     U, S, V = @constinferred svd_compact(A)
     @test U * S * V == A
@@ -220,6 +196,22 @@ end
     # @test Q === A
     # @test R === A
 
+    # A = Eye(4, 3)
+    # Q, R = @constinferred left_polar(A)
+    # @test Q * R == A
+    # @test size(Q) == (4, 3)
+    # @test size(R) == (3, 3)
+    # @test Q == Matrix(I, (4, 3))
+    # @test Q isa Eye
+    # @test R == I
+    # @test R isa Eye
+
+    # A = Eye(3)
+    # Q, R = @constinferred left_polar(A)
+    # @test Q * R == A
+    # @test Q === A
+    # @test R === A
+
     A = Eye(3, 4)
     L, Q = @constinferred lq_compact(A)
     @test L * Q == A
@@ -252,6 +244,22 @@ end
     @test L === A
     @test Q === A
 
+    # A = Eye(3, 4)
+    # L, Q = @constinferred right_polar(A)
+    # @test L * Q == A
+    # @test size(L) == (3, 3)
+    # @test size(Q) == (3, 4)
+    # @test L == I
+    # @test L isa Eye
+    # @test Q == Matrix(I, (3, 4))
+    # @test Q isa Eye
+
+    # A = Eye(3)
+    # L, Q = @constinferred right_polar(A)
+    # @test L * Q == A
+    # @test L === A
+    # @test Q === A
+
     A = Eye(3, 4)
     U, S, V = @constinferred svd_compact(A)
     @test U * S * V == A

From 880e8d5a9e5be66e1b1915576c5c440b01cc1295 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 2 Jun 2025 08:46:13 -0400
Subject: [PATCH 08/13] Trunc Eye

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 42 +++++++++++++++++++++++++++-
 test/fillarrays.jl                   | 33 ++++++++++++++++++++--
 2 files changed, 72 insertions(+), 3 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 62f729bf..9d281d38 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -2,7 +2,9 @@ module MatrixAlgebraKitFillArraysExt
 
 using LinearAlgebra
 using MatrixAlgebraKit
-using MatrixAlgebraKit: AbstractAlgorithm, check_input, diagview
+using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, TruncationStrategy,
+                        check_input, diagview,
+                        findtruncated, select_algorithm, select_truncation
 using FillArrays
 using FillArrays: AbstractZerosMatrix, OnesVector, RectDiagonal, SquareEye
 
@@ -177,6 +179,26 @@ for f in [:eig_full!, :eigh_full!]
     end
 end
 
+for f in [:eig_trunc!, :eigh_trunc!]
+    @eval begin
+        # TODO: Delete this when `select_algorithm` is generalized.
+        function MatrixAlgebraKit.select_algorithm(::typeof($f), ::Type{A}, alg;
+                                                   trunc=nothing,
+                                                   kwargs...) where {A<:Eye}
+            alg_eig = select_algorithm(eig_full!, A, alg; kwargs...)
+            return TruncatedAlgorithm(alg_eig, select_truncation(trunc))
+        end
+        # TODO: I think it would be better to dispatch on the algorithm here,
+        # rather than the output types.
+        function MatrixAlgebraKit.truncate!(::typeof($f), (D, V)::Tuple{Eye,Eye},
+                                            strategy::TruncationStrategy)
+            ind = findtruncated(diagview(D), strategy)
+            return Diagonal(diagview(D)[ind]),
+                   Eye((axes(V, 1), only(axes(axes(V, 2)[ind]))))
+        end
+    end
+end
+
 for f in [:eig_vals!, :eigh_vals!]
     @eval begin
         function MatrixAlgebraKit.check_input(::typeof($f), A::Eye, F)
@@ -262,6 +284,24 @@ function MatrixAlgebraKit.svd_full!(A::SquareEye, F, alg::EyeAlgorithm)
     return (A, A, A)
 end
 
+# TODO: Delete this when `select_algorithm` is generalized.
+function MatrixAlgebraKit.select_algorithm(::typeof(svd_trunc!), ::Type{A}, alg;
+                                           trunc=nothing,
+                                           kwargs...) where {A<:Eye}
+    alg_eig = select_algorithm(eig_full!, A, alg; kwargs...)
+    return TruncatedAlgorithm(alg_eig, select_truncation(trunc))
+end
+# TODO: I think it would be better to dispatch on the algorithm here,
+# rather than the output types.
+function MatrixAlgebraKit.truncate!(::typeof(svd_trunc!), (U, S, V)::Tuple{Eye,Eye,Eye},
+                                    strategy::TruncationStrategy)
+    ind = findtruncated(diagview(S), strategy)
+    U′ = Eye((axes(U, 1), only(axes(axes(U, 2)[ind]))))
+    S′ = Diagonal(diagview(S)[ind])
+    V′ = Eye((only(axes(axes(V, 1)[ind])), axes(V, 2)))
+    return (U′, S′, V′)
+end
+
 function MatrixAlgebraKit.svd_vals!(A::Eye, F, alg::EyeAlgorithm)
     return diagview(A)
 end
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index f772d641..1baf9fe6 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -3,13 +3,14 @@ using LinearAlgebra
 using Test
 using TestExtras
 using FillArrays
+using FillArrays: SquareEye
 
 @testset "Zeros" begin
     for f in [:eig_full, :eigh_full]
         @eval begin
             A = Zeros(3, 3)
             D, V = @constinferred $f(A)
-            @test A * V == D * V
+            @test A * V == V * D
             @test size(D) == size(A)
             @test size(V) == size(A)
             @test iszero(D)
@@ -142,7 +143,7 @@ end
         @eval begin
             for A in (Eye(3), Eye(3, 3))
                 local D, V = @constinferred $f(A)
-                @test A * V == D * V
+                @test A * V == V * D
                 @test size(D) == size(A)
                 @test size(V) == size(A)
                 @test V == I
@@ -153,6 +154,21 @@ end
         end
     end
 
+    for f in [:eig_trunc, :eigh_trunc]
+        @eval begin
+            for A in (Eye(3), Eye(3, 3))
+                local D, V = @constinferred $f(A; trunc=(; maxrank=2))
+                @test A * V == V * D
+                @test size(D) == (2, 2)
+                @test size(V) == (3, 2)
+                @test D == Eye(2, 2)
+                @test D isa SquareEye
+                @test V == Eye(3, 2)
+                @test V isa Eye
+            end
+        end
+    end
+
     for f in [:eig_vals, :eigh_vals]
         @eval begin
             for A in (Eye(3), Eye(3, 3))
@@ -300,6 +316,19 @@ end
     @test S === A
     @test V === A
 
+    A = Eye(3, 4)
+    U, S, V = @constinferred svd_trunc(A; trunc=(; maxrank=2))
+    @test U * S * V == Eye(3, 2) * Eye(2, 2) * Eye(2, 4)
+    @test size(U) == (3, 2)
+    @test size(S) == (2, 2)
+    @test size(V) == (2, 4)
+    @test S == Eye(2, 2)
+    @test S isa Eye
+    @test U == Eye(3, 2)
+    @test U isa Eye
+    @test V == Eye(2, 4)
+    @test V isa Eye
+
     A = Eye(3, 4)
     D = @constinferred svd_vals(A)
     @test size(D) == (minimum(size(A)),)

From ed8d2bd28c633e7300c638ac71abc4ee01ed7e7a Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 2 Jun 2025 09:43:11 -0400
Subject: [PATCH 09/13] More support for truncation

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 44 ++++++++++++++++++++++++++--
 test/fillarrays.jl                   | 27 +++++++++++++++++
 2 files changed, 69 insertions(+), 2 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 9d281d38..e6209173 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -60,6 +60,27 @@ for f in [:eig_full!, :eigh_full!]
     end
 end
 
+for f in [:eig_trunc!, :eigh_trunc!]
+    @eval begin
+        # TODO: Delete this when `select_algorithm` is generalized.
+        function MatrixAlgebraKit.select_algorithm(::typeof($f), ::Type{A}, alg;
+                                                   trunc=nothing,
+                                                   kwargs...) where {A<:Zeros}
+            alg_eig = select_algorithm(eig_full!, A, alg; kwargs...)
+            return TruncatedAlgorithm(alg_eig, select_truncation(trunc))
+        end
+        # TODO: I think it would be better to dispatch on the algorithm here,
+        # rather than the output types.
+        function MatrixAlgebraKit.truncate!(::typeof($f), (D, V)::Tuple{Zeros,Eye},
+                                            strategy::TruncationStrategy)
+            ind = findtruncated(diagview(D), strategy)
+            D′ = D[ind, ind]
+            V′ = Eye((axes(V, 1), only(axes(axes(V, 2)[ind]))))
+            return (D′, V′)
+        end
+    end
+end
+
 for f in [:eig_vals!, :eigh_vals!]
     @eval begin
         function MatrixAlgebraKit.check_input(::typeof($f), A::AbstractZerosMatrix, F)
@@ -127,6 +148,24 @@ function MatrixAlgebraKit.svd_full!(A::AbstractZerosMatrix, F, alg::ZerosAlgorit
     return (Eye((ax[1], ax[1])), Zeros(ax), Eye((ax[2], ax[2])))
 end
 
+# TODO: Delete this when `select_algorithm` is generalized.
+function MatrixAlgebraKit.select_algorithm(::typeof(svd_trunc!), ::Type{A}, alg;
+                                           trunc=nothing,
+                                           kwargs...) where {A<:Zeros}
+    alg_eig = select_algorithm(eig_full!, A, alg; kwargs...)
+    return TruncatedAlgorithm(alg_eig, select_truncation(trunc))
+end
+# TODO: I think it would be better to dispatch on the algorithm here,
+# rather than the output types.
+function MatrixAlgebraKit.truncate!(::typeof(svd_trunc!), (U, S, V)::Tuple{Eye,Zeros,Eye},
+                                    strategy::TruncationStrategy)
+    ind = findtruncated(diagview(S), strategy)
+    U′ = Eye((axes(U, 1), only(axes(axes(U, 2)[ind]))))
+    S′ = S[ind, ind]
+    V′ = Eye((only(axes(axes(V, 1)[ind])), axes(V, 2)))
+    return (U′, S′, V′)
+end
+
 function MatrixAlgebraKit.svd_vals!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
     return diagview(A)
 end
@@ -193,8 +232,9 @@ for f in [:eig_trunc!, :eigh_trunc!]
         function MatrixAlgebraKit.truncate!(::typeof($f), (D, V)::Tuple{Eye,Eye},
                                             strategy::TruncationStrategy)
             ind = findtruncated(diagview(D), strategy)
-            return Diagonal(diagview(D)[ind]),
-                   Eye((axes(V, 1), only(axes(axes(V, 2)[ind]))))
+            D′ = Diagonal(diagview(D)[ind])
+            U′ = Eye((axes(V, 1), only(axes(axes(V, 2)[ind]))))
+            return (D′, U′)
         end
     end
 end
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index 1baf9fe6..16883c48 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -20,6 +20,20 @@ using FillArrays: SquareEye
         end
     end
 
+    for f in [:eig_trunc, :eigh_trunc]
+        @eval begin
+            A = Zeros(3, 3)
+            D, V = @constinferred $f(A; trunc=(; maxrank=2))
+            @test A * V == V * D
+            @test size(D) == (2, 2)
+            @test size(V) == (3, 2)
+            @test D == Zeros(2, 2)
+            @test D isa Zeros
+            @test V == Eye(3, 2)
+            @test V isa Eye
+        end
+    end
+
     for f in [:eig_vals, :eigh_vals]
         @eval begin
             A = Zeros(3, 3)
@@ -131,6 +145,19 @@ using FillArrays: SquareEye
     @test V == I
     @test V isa Eye
 
+    A = Zeros(3, 4)
+    U, S, V = @constinferred svd_trunc(A; trunc=(; maxrank=2))
+    @test U * S * V == Eye(3, 2) * Zeros(2, 2) * Eye(2, 4)
+    @test size(U) == (3, 2)
+    @test size(S) == (2, 2)
+    @test size(V) == (2, 4)
+    @test S == Zeros(2, 2)
+    @test S isa Zeros
+    @test U == Eye(3, 2)
+    @test U isa Eye
+    @test V == Eye(2, 4)
+    @test V isa Eye
+
     A = Zeros(3, 4)
     D = @constinferred svd_vals(A)
     @test size(D) == (minimum(size(A)),)

From fa786cbe4f0a674c6d99ed1553b2efaadd00c618 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 4 Jun 2025 19:03:13 -0400
Subject: [PATCH 10/13] Reenable fixed QR tests

---
 Project.toml       |  2 +-
 test/fillarrays.jl | 36 ++++++++++++++++++------------------
 2 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/Project.toml b/Project.toml
index 28c8db66..aa30c584 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MatrixAlgebraKit"
 uuid = "6c742aac-3347-4629-af66-fc926824e5e4"
 authors = ["Jutho <jutho.haegeman@ugent.be> and contributors"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index 16883c48..e9edb52f 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -44,25 +44,25 @@ using FillArrays: SquareEye
         end
     end
 
-    # A = Zeros(4, 3)
-    # Q, R = @constinferred qr_compact(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 3)
-    # @test size(R) == (3, 3)
-    # @test Q == Matrix(I, (4, 3))
-    # @test Q isa Eye
-    # @test iszero(R)
-    # @test R isa Zeros
+    A = Zeros(4, 3)
+    Q, R = @constinferred qr_compact(A)
+    @test Q * R == A
+    @test size(Q) == (4, 3)
+    @test size(R) == (3, 3)
+    @test Q == Matrix(I, (4, 3))
+    @test Q isa Eye
+    @test iszero(R)
+    @test R isa Zeros
 
-    # A = Zeros(4, 3)
-    # Q, R = @constinferred qr_full(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 4)
-    # @test size(R) == (4, 3)
-    # @test Q == I
-    # @test Q isa Eye
-    # @test iszero(R)
-    # @test R isa Zeros
+    A = Zeros(4, 3)
+    Q, R = @constinferred qr_full(A)
+    @test Q * R == A
+    @test size(Q) == (4, 4)
+    @test size(R) == (4, 3)
+    @test Q == I
+    @test Q isa Eye
+    @test iszero(R)
+    @test R isa Zeros
 
     # A = Zeros(4, 3)
     # Q, R = @constinferred left_polar(A)

From c846fa0544adbf7e9096333f1d3b6579190be748 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 4 Jun 2025 19:22:47 -0400
Subject: [PATCH 11/13] Fix polar

---
 ext/MatrixAlgebraKitFillArraysExt.jl | 26 +++++++++++++++++++-
 test/fillarrays.jl                   | 36 ++++++++++++++--------------
 2 files changed, 43 insertions(+), 19 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index e6209173..7c1c34a7 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -14,7 +14,7 @@ end
 
 struct ZerosAlgorithm <: AbstractAlgorithm end
 
-for f in [:eig, :eigh, :lq, :qr, :svd]
+for f in [:eig, :eigh, :lq, :polar, :qr, :svd]
     ff = Symbol("default_", f, "_algorithm")
     @eval begin
         function MatrixAlgebraKit.$ff(::Type{<:AbstractZerosMatrix}; kwargs...)
@@ -170,6 +170,30 @@ function MatrixAlgebraKit.svd_vals!(A::AbstractZerosMatrix, F, alg::ZerosAlgorit
     return diagview(A)
 end
 
+function MatrixAlgebraKit.check_input(::typeof(left_polar!), A::AbstractZerosMatrix, F)
+    m, n = size(A)
+    m >= n ||
+        throw(ArgumentError("input matrix needs at least as many rows as columns"))
+    return nothing
+end
+function MatrixAlgebraKit.left_polar!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    check_input(left_polar!, A, F)
+    U, S, Vᴴ = svd_compact(A)
+    return (Eye((axes(U, 1), axes(Vᴴ, 2))), Vᴴ' * S * Vᴴ)
+end
+
+function MatrixAlgebraKit.check_input(::typeof(right_polar!), A::AbstractZerosMatrix, F)
+    m, n = size(A)
+    n >= m ||
+        throw(ArgumentError("input matrix needs at least as many columns as rows"))
+    return nothing
+end
+function MatrixAlgebraKit.right_polar!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    check_input(right_polar!, A, F)
+    U, S, Vᴴ = svd_compact(A)
+    return (U * S * U', Eye((axes(U, 1), axes(Vᴴ, 2))))
+end
+
 struct EyeAlgorithm <: AbstractAlgorithm end
 
 for f in [:eig, :eigh, :lq, :qr, :polar, :svd]
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index e9edb52f..a69faa8d 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -64,15 +64,15 @@ using FillArrays: SquareEye
     @test iszero(R)
     @test R isa Zeros
 
-    # A = Zeros(4, 3)
-    # Q, R = @constinferred left_polar(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 3)
-    # @test size(R) == (3, 3)
-    # @test Q == Matrix(I, (4, 3))
-    # @test Q isa Eye
-    # @test iszero(R)
-    # @test R isa Zeros
+    A = Zeros(4, 3)
+    Q, R = @constinferred left_polar(A)
+    @test Q * R == A
+    @test size(Q) == (4, 3)
+    @test size(R) == (3, 3)
+    @test Q == Matrix(I, (4, 3))
+    @test Q isa Eye
+    @test iszero(R)
+    @test R isa Zeros
 
     A = Zeros(3, 4)
     L, Q = @constinferred lq_compact(A)
@@ -109,15 +109,15 @@ using FillArrays: SquareEye
     @test Q == Eye(4)
     @test Q isa Eye
 
-    # A = Zeros(3, 4)
-    # L, Q = @constinferred right_polar(A)
-    # @test L * Q == A
-    # @test size(L) == (3, 3)
-    # @test size(Q) == (3, 4)
-    # @test iszero(L)
-    # @test L isa Zeros
-    # @test Q == Matrix(I, (3, 4))
-    # @test Q isa Eye
+    A = Zeros(3, 4)
+    L, Q = @constinferred right_polar(A)
+    @test L * Q == A
+    @test size(L) == (3, 3)
+    @test size(Q) == (3, 4)
+    @test iszero(L)
+    @test L isa Zeros
+    @test Q == Matrix(I, (3, 4))
+    @test Q isa Eye
 
     A = Zeros(3, 4)
     U, S, V = @constinferred svd_compact(A)

From f16ebfc193f0ecdaacac04554b6e575afdf67172 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 4 Jun 2025 19:24:50 -0400
Subject: [PATCH 12/13] Test left_orth, right_orth

---
 test/fillarrays.jl | 36 ++++++++++++++++++++----------------
 1 file changed, 20 insertions(+), 16 deletions(-)

diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index a69faa8d..c2879c22 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -44,15 +44,17 @@ using FillArrays: SquareEye
         end
     end
 
-    A = Zeros(4, 3)
-    Q, R = @constinferred qr_compact(A)
-    @test Q * R == A
-    @test size(Q) == (4, 3)
-    @test size(R) == (3, 3)
-    @test Q == Matrix(I, (4, 3))
-    @test Q isa Eye
-    @test iszero(R)
-    @test R isa Zeros
+    for f in (qr_compact, left_orth)
+        A = Zeros(4, 3)
+        Q, R = @constinferred f(A)
+        @test Q * R == A
+        @test size(Q) == (4, 3)
+        @test size(R) == (3, 3)
+        @test Q == Matrix(I, (4, 3))
+        @test Q isa Eye
+        @test iszero(R)
+        @test R isa Zeros
+    end
 
     A = Zeros(4, 3)
     Q, R = @constinferred qr_full(A)
@@ -84,13 +86,15 @@ using FillArrays: SquareEye
     @test Q == Matrix(I, (3, 4))
     @test Q isa Eye
 
-    A = Zeros(3, 4)
-    L, Q = @constinferred lq_compact(A)
-    @test L * Q == A
-    @test L == Zeros(3, 3)
-    @test L isa Zeros
-    @test Q == Eye(3, 4)
-    @test Q isa Eye
+    for f in (lq_compact, right_orth)
+        A = Zeros(3, 4)
+        L, Q = @constinferred lq_compact(A)
+        @test L * Q == A
+        @test L == Zeros(3, 3)
+        @test L isa Zeros
+        @test Q == Eye(3, 4)
+        @test Q isa Eye
+    end
 
     A = Zeros(3, 4)
     L, Q = @constinferred lq_full(A)

From 34200fd70d851915b313ef5f16dcb0a715914a21 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 4 Jun 2025 22:07:35 -0400
Subject: [PATCH 13/13] left and right polar for Eye

---
 ext/MatrixAlgebraKitFillArraysExt.jl |  30 ++++-
 test/fillarrays.jl                   | 162 +++++++++++++--------------
 2 files changed, 107 insertions(+), 85 deletions(-)

diff --git a/ext/MatrixAlgebraKitFillArraysExt.jl b/ext/MatrixAlgebraKitFillArraysExt.jl
index 7c1c34a7..8d09c012 100644
--- a/ext/MatrixAlgebraKitFillArraysExt.jl
+++ b/ext/MatrixAlgebraKitFillArraysExt.jl
@@ -179,7 +179,7 @@ end
 function MatrixAlgebraKit.left_polar!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
     check_input(left_polar!, A, F)
     U, S, Vᴴ = svd_compact(A)
-    return (Eye((axes(U, 1), axes(Vᴴ, 2))), Vᴴ' * S * Vᴴ)
+    return (Eye((axes(A, 1), axes(A, 2))), Vᴴ' * S * Vᴴ)
 end
 
 function MatrixAlgebraKit.check_input(::typeof(right_polar!), A::AbstractZerosMatrix, F)
@@ -370,4 +370,32 @@ function MatrixAlgebraKit.svd_vals!(A::Eye, F, alg::EyeAlgorithm)
     return diagview(A)
 end
 
+function MatrixAlgebraKit.check_input(::typeof(left_polar!), A::Eye, F)
+    m, n = size(A)
+    m >= n ||
+        throw(ArgumentError("input matrix needs at least as many rows as columns"))
+    return nothing
+end
+function MatrixAlgebraKit.left_polar!(A::Eye, F, alg::EyeAlgorithm)
+    check_input(left_polar!, A, F)
+    return (Eye((axes(A, 1), axes(A, 2))), Eye((axes(A, 2), axes(A, 2))))
+end
+function MatrixAlgebraKit.left_polar!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A)
+end
+
+function MatrixAlgebraKit.check_input(::typeof(right_polar!), A::Eye, F)
+    m, n = size(A)
+    n >= m ||
+        throw(ArgumentError("input matrix needs at least as many columns as rows"))
+    return nothing
+end
+function MatrixAlgebraKit.right_polar!(A::Eye, F, alg::EyeAlgorithm)
+    check_input(right_polar!, A, F)
+    return (Eye((axes(A, 1), axes(A, 1))), Eye((axes(A, 1), axes(A, 2))))
+end
+function MatrixAlgebraKit.right_polar!(A::SquareEye, F, alg::EyeAlgorithm)
+    return (A, A)
+end
+
 end
diff --git a/test/fillarrays.jl b/test/fillarrays.jl
index c2879c22..ed278330 100644
--- a/test/fillarrays.jl
+++ b/test/fillarrays.jl
@@ -76,19 +76,9 @@ using FillArrays: SquareEye
     @test iszero(R)
     @test R isa Zeros
 
-    A = Zeros(3, 4)
-    L, Q = @constinferred lq_compact(A)
-    @test L * Q == A
-    @test size(L) == (3, 3)
-    @test size(Q) == (3, 4)
-    @test iszero(L)
-    @test L isa Zeros
-    @test Q == Matrix(I, (3, 4))
-    @test Q isa Eye
-
     for f in (lq_compact, right_orth)
         A = Zeros(3, 4)
-        L, Q = @constinferred lq_compact(A)
+        L, Q = @constinferred f(A)
         @test L * Q == A
         @test L == Zeros(3, 3)
         @test L isa Zeros
@@ -211,69 +201,73 @@ end
         end
     end
 
-    # A = Eye(4, 3)
-    # Q, R = @constinferred qr_compact(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 3)
-    # @test size(R) == (3, 3)
-    # @test Q == Matrix(I, (4, 3))
-    # @test Q isa Eye
-    # @test R == I
-    # @test R isa Eye
-
-    # A = Eye(3)
-    # Q, R = @constinferred qr_compact(A)
-    # @test Q * R == A
-    # @test Q === A
-    # @test R === A
-
-    # A = Eye(4, 3)
-    # Q, R = @constinferred qr_full(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 4)
-    # @test size(R) == (4, 3)
-    # @test Q == I
-    # @test Q isa Eye
-    # @test R == I
-    # @test R isa Eye
-
-    # A = Eye(3)
-    # Q, R = @constinferred qr_full(A)
-    # @test Q * R == A
-    # @test Q === A
-    # @test R === A
-
-    # A = Eye(4, 3)
-    # Q, R = @constinferred left_polar(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 3)
-    # @test size(R) == (3, 3)
-    # @test Q == Matrix(I, (4, 3))
-    # @test Q isa Eye
-    # @test R == I
-    # @test R isa Eye
-
-    # A = Eye(3)
-    # Q, R = @constinferred left_polar(A)
-    # @test Q * R == A
-    # @test Q === A
-    # @test R === A
+    for f in (qr_compact, left_orth)
+        A = Eye(4, 3)
+        Q, R = @constinferred f(A)
+        @test Q * R == A
+        @test size(Q) == (4, 3)
+        @test size(R) == (3, 3)
+        @test Q == Matrix(I, (4, 3))
+        @test Q isa Eye
+        @test R == I
+        @test R isa Eye
 
-    A = Eye(3, 4)
-    L, Q = @constinferred lq_compact(A)
-    @test L * Q == A
-    @test size(L) == (3, 3)
-    @test size(Q) == (3, 4)
-    @test L == I
-    @test L isa Eye
-    @test Q == Matrix(I, (3, 4))
+        A = Eye(3)
+        Q, R = @constinferred f(A)
+        @test Q * R == A
+        @test Q === A
+        @test R === A
+    end
+
+    A = Eye(4, 3)
+    Q, R = @constinferred qr_full(A)
+    @test Q * R == A
+    @test size(Q) == (4, 4)
+    @test size(R) == (4, 3)
+    @test Q == I
     @test Q isa Eye
+    @test R == Eye(4, 3)
+    @test R isa Eye
 
     A = Eye(3)
-    L, Q = @constinferred lq_compact(A)
-    @test L * Q == A
-    @test L === A
+    Q, R = @constinferred qr_full(A)
+    @test Q * R == A
     @test Q === A
+    @test R === A
+
+    A = Eye(4, 3)
+    Q, R = @constinferred left_polar(A)
+    @test Q * R == A
+    @test size(Q) == (4, 3)
+    @test size(R) == (3, 3)
+    @test Q == Matrix(I, (4, 3))
+    @test Q isa Eye
+    @test R == I
+    @test R isa Eye
+
+    A = Eye(3)
+    Q, R = @constinferred left_polar(A)
+    @test Q * R == A
+    @test Q === A
+    @test R === A
+
+    for f in (lq_compact, right_orth)
+        A = Eye(3, 4)
+        L, Q = @constinferred lq_compact(A)
+        @test L * Q == A
+        @test size(L) == (3, 3)
+        @test size(Q) == (3, 4)
+        @test L == I
+        @test L isa Eye
+        @test Q == Matrix(I, (3, 4))
+        @test Q isa Eye
+
+        A = Eye(3)
+        L, Q = @constinferred lq_compact(A)
+        @test L * Q == A
+        @test L === A
+        @test Q === A
+    end
 
     A = Eye(3, 4)
     L, Q = @constinferred lq_full(A)
@@ -291,21 +285,21 @@ end
     @test L === A
     @test Q === A
 
-    # A = Eye(3, 4)
-    # L, Q = @constinferred right_polar(A)
-    # @test L * Q == A
-    # @test size(L) == (3, 3)
-    # @test size(Q) == (3, 4)
-    # @test L == I
-    # @test L isa Eye
-    # @test Q == Matrix(I, (3, 4))
-    # @test Q isa Eye
-
-    # A = Eye(3)
-    # L, Q = @constinferred right_polar(A)
-    # @test L * Q == A
-    # @test L === A
-    # @test Q === A
+    A = Eye(3, 4)
+    L, Q = @constinferred right_polar(A)
+    @test L * Q == A
+    @test size(L) == (3, 3)
+    @test size(Q) == (3, 4)
+    @test L == I
+    @test L isa Eye
+    @test Q == Matrix(I, (3, 4))
+    @test Q isa Eye
+
+    A = Eye(3)
+    L, Q = @constinferred right_polar(A)
+    @test L * Q == A
+    @test L === A
+    @test Q === A
 
     A = Eye(3, 4)
     U, S, V = @constinferred svd_compact(A)