Write isapprox in terms of isapprox of the factors

Author: mtfishman

src/kroneckerarray.jl

@@ -378,26 +378,29 @@ end
 
 using LinearAlgebra: dot, promote_leaf_eltypes
 function Base.isapprox(
-        a::AbstractKroneckerArray, b::AbstractKroneckerArray;
-        atol::Real = 0,
+        a::AbstractKroneckerArray, b::AbstractKroneckerArray; atol::Real = 0,
         rtol::Real = Base.rtoldefault(promote_leaf_eltypes(a), promote_leaf_eltypes(b), atol),
-        norm::Function = norm
     )
     a1, a2 = arg1(a), arg2(a)
     b1, b2 = arg1(b), arg2(b)
-    d = if a1 == b1
-        norm(a1) * norm(a2 - b2)
+    if a1 == b1
+        return isapprox(a2, b2; atol = atol / norm(a1), rtol)
     elseif a2 == b2
-        norm(a1 - b1) * norm(b2)
+        return isapprox(a1, b1; atol = atol / norm(a2), rtol)
     else
-        # This could be defined as `KroneckerArrays.dist_kronecker(a, b)`, but that might have
-        # numerical precision issues so for now we just error.
-        error(
-            "`isapprox` not implemented for KroneckerArrays where both arguments differ. " *
-                "In those cases, you can use `isapprox(collect(a), collect(b); kwargs...)`."
+        # This could be defined as:
+        # ```julia
+        # d = KroneckerArrays.dist_kronecker(a, b)
+        # iszero(rtol) ? d <= atol : d <= max(atol, rtol * max(norm(a), norm(b)))
+        # ```
+        # but that might have numerical precision issues so for now we just error.
+        throw(
+            ArgumentError(
+                "`isapprox` not implemented for KroneckerArrays where both arguments differ. " *
+                    "In those cases, you can use `isapprox(collect(a), collect(b); kwargs...)`."
+            )
         )
     end
-    return iszero(rtol) ? d <= atol : d <= max(atol, rtol * max(norm(a), norm(b)))
 end
 
 function Base.iszero(a::AbstractKroneckerArray)

src/linearalgebra.jl

@@ -58,7 +58,7 @@ function LinearAlgebra.tr(a::AbstractKroneckerArray)
 end
 
 using LinearAlgebra: norm
-function LinearAlgebra.norm(a::AbstractKroneckerArray, p::Int = 2)
+function LinearAlgebra.norm(a::AbstractKroneckerArray, p::Real = 2)
     return norm(arg1(a), p) * norm(arg2(a), p)
 end
 

test/test_basics.jl

@@ -237,9 +237,47 @@ elts = (Float32, Float64, ComplexF32, ComplexF64)
         end
     end
 
-    # KroneckerArrays.dist
+    # isapprox
+
+    rng = StableRNG(123)
+    a1 = randn(rng, elt, (2, 2))
+    a = a1 ⊗ randn(rng, elt, (3, 3))
+    b = a1 ⊗ randn(rng, elt, (3, 3))
+    @test isapprox(a, b; atol = norm(a - b) * (1 + 2eps(real(elt))))
+    @test !isapprox(a, b; atol = norm(a - b) * (1 - 2eps(real(elt))))
+    @test isapprox(
+        a, b;
+        rtol = norm(a - b) / max(norm(a), norm(b)) * (1 + 2eps(real(elt)))
+    )
+    @test !isapprox(
+        a, b;
+        rtol = norm(a - b) / max(norm(a), norm(b)) * (1 - 2eps(real(elt)))
+    )
+    @test isapprox(
+        a, b; atol = norm(a - b) * (1 + 2eps(real(elt))),
+        rtol = norm(a - b) / max(norm(a), norm(b)) * (1 + 2eps(real(elt)))
+    )
+    @test isapprox(
+        a, b; atol = norm(a - b) * (1 + 2eps(real(elt))),
+        rtol = norm(a - b) / max(norm(a), norm(b)) * (1 - 2eps(real(elt)))
+    )
+    @test isapprox(
+        a, b; atol = norm(a - b) * (1 - 2eps(real(elt))),
+        rtol = norm(a - b) / max(norm(a), norm(b)) * (1 + 2eps(real(elt)))
+    )
+    @test !isapprox(
+        a, b; atol = norm(a - b) * (1 - 2eps(real(elt))),
+        rtol = norm(a - b) / max(norm(a), norm(b)) * (1 - 2eps(real(elt)))
+    )
+
+    a = randn(elt, (2, 2)) ⊗ randn(elt, (3, 3))
+    b = randn(elt, (2, 2)) ⊗ randn(elt, (3, 3))
+    @test_throws ArgumentError isapprox(a, b)
+
+    # KroneckerArrays.dist_kronecker
     rng = StableRNG(123)
-    a = randn(rng, 100, 100) ⊗ randn(rng, 100, 100)
-    b = (arg1(a) + 1.0e-1 * randn(rng, size(arg1(a)))) ⊗ (arg2(a) + 1.0e-1 * randn(rng, size(arg2(a))))
+    a = randn(rng, (100, 100)) ⊗ randn(rng, (100, 100))
+    b = (arg1(a) + randn(rng, size(arg1(a))) / 10) ⊗
+        (arg2(a) + randn(rng, size(arg2(a))) / 10)
     @test KroneckerArrays.dist_kronecker(a, b) ≈ norm(collect(a) - collect(b)) rtol = 1.0e-2
 end