Merge pull request #45 from SciML/zeromatrix

add zeromatrix

Author: ChrisRackauckas

Project.toml

@@ -1,6 +1,6 @@
 name = "ArrayInterface"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "2.6.2"
+version = "2.7.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

README.md

@@ -69,6 +69,22 @@ cheaply.
 
 Return an instance of the LU factorization object with the correct type cheaply.
 
+## safevec(v)
+
+Is a form of `vec` which is safe for all values in vector spaces, i.e. if
+is already a vector, like an AbstractVector or Number, it will return said
+AbstractVector or Number.
+
+## zeromatrix(u::AbstractVector)
+
+Creates the zero'd matrix version of `u`. Note that this is unique because
+`similar(u,length(u),length(u))` returns a mutable type, so is not type-matching,
+while `fill(zero(eltype(u)),length(u),length(u))` doesn't match the array type,
+i.e. you'll get a CPU array from a GPU array. The generic fallback is
+`u .* u' .* false` which works on a surprising number of types, but can be broken
+with weird (recursive) broadcast overloads. For higher order tensors, this
+returns the matrix linear operator type which acts on the `vec` of the array.
+
 ## List of things to add
 
 - https://github.com/JuliaLang/julia/issues/22216

src/ArrayInterface.jl

@@ -420,6 +420,33 @@ Return the number.
 """
 lu_instance(a::Number) = a
 
+"""
+safevec(v)
+
+Is a form of `vec` which is safe for all values in vector spaces, i.e. if
+is already a vector, like an AbstractVector or Number, it will return said
+AbstractVector or Number.
+"""
+safevec(v) = vec(v)
+safevec(v::Number) = v
+safevec(v::AbstractVector) = v
+
+"""
+zeromatrix(u::AbstractVector)
+
+Creates the zero'd matrix version of `u`. Note that this is unique because
+`similar(u,length(u),length(u))` returns a mutable type, so is not type-matching,
+while `fill(zero(eltype(u)),length(u),length(u))` doesn't match the array type,
+i.e. you'll get a CPU array from a GPU array. The generic fallback is
+`u .* u' .* false` which works on a surprising number of types, but can be broken
+with weird (recursive) broadcast overloads. For higher order tensors, this
+returns the matrix linear operator type which acts on the `vec` of the array.
+"""
+function zeromatrix(u)
+  x = safevec(u)
+  x .* x' .* false
+end
+
 function __init__()
 
   @require SuiteSparse="4607b0f0-06f3-5cda-b6b1-a6196a1729e9" begin

test/runtests.jl

@@ -159,3 +159,6 @@ using ArrayInterface: issingular
         @test all(!issingular, [UnitLowerTriangular(R), UnitUpperTriangular(R), UnitUpperTriangular(R)'])
     end
 end
+
+using ArrayInterface: zeromatrix
+@test zeromatrix(rand(4,4,4)) == zeros(4*4*4,4*4*4)