Merge pull request #113 from Lewih/ArrayPartition-New-LA-overloads

Added linalg mul! overloads for ArrayPartition with tests

Author: ChrisRackauckas

src/array_partition.jl

@@ -370,3 +370,29 @@ function LinearAlgebra._swap_rows!(B::ArrayPartition, i::Integer, j::Integer)
     B[i], B[j] = B[j], B[i]
     return B
 end
+
+# linalg mul! overloads for ArrayPartition
+function LinearAlgebra.mul!(C::ArrayPartition, A::ArrayPartition, B::AbstractArray)
+    if length(C.x) != length(A.x)
+        throw(DimensionMismatch("Length of C, $(length(C.x)), does not match length of A, $(length(A.x))"))
+    end
+
+    for index = 1:length(C.x)
+        mul!(C.x[index], A.x[index], B)
+    end
+    return C
+end
+
+function LinearAlgebra.mul!(C::ArrayPartition, A::ArrayPartition, B::ArrayPartition)
+    if length(C.x) != length(A.x)
+        throw(DimensionMismatch("Length of C, $(length(C.x)), does not match length of A, $(length(B.x))"))
+    end
+    if length(A.x) != length(B.x)
+        throw(DimensionMismatch("Length of A, $(length(A.x)), does not match length of B, $(length(B.x))"))
+    end
+
+    for index = 1:length(C.x)
+        mul!(C.x[index], A.x[index], B.x[index])
+    end
+    return C
+end

test/linalg.jl

@@ -8,7 +8,7 @@ b = ArrayPartition(bb)
 A = randn(MersenneTwister(123), n+m, n+m)
 
 for T in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
-    B = T(A)
+    local B = T(A)
     @test B*Array(B \ b) ≈ b
     bbb = copy(b)
     @test ldiv!(bbb, B, b) === bbb
@@ -26,3 +26,24 @@ for ff in (lu, svd, qr)
     @test ldiv!(FF, bbb) === bbb
     @test A*bbb ≈ b
 end
+
+# linalg mul! overloads
+n, m, l = 5, 6, 7
+bb = rand(n, n), rand(m, n), rand(l, n)
+cc = rand(n), rand(n), rand(n)
+dd = rand(n), rand(m), rand(l)
+b = ArrayPartition(bb)
+c = ArrayPartition(cc)
+d = ArrayPartition(dd)
+A = rand(n)
+for T in (Array{Float64}, Array{ComplexF64},)
+    local B = T(A)
+    mul!(d, b, A)
+    for i = 1:length(c.x)
+        @test d.x[i] == b.x[i] * A
+    end
+    mul!(d, b, c)
+    for i = 1:length(d.x)
+        @test d.x[i] == b.x[i] * c.x[i]
+    end
+end

test/partitions_test.jl

@@ -111,8 +111,8 @@ S = [
     ]
 
 for sizes in S
-  x = ArrayPartition( randn.(sizes[1]) )
-  y = ArrayPartition( zeros.(sizes[2]) )
+  local x = ArrayPartition( randn.(sizes[1]) )
+  local y = ArrayPartition( zeros.(sizes[2]) )
   y_array = zeros(length(x))
   copyto!(y,x)           #testing Base.copyto!(dest::ArrayPartition,A::ArrayPartition)
   copyto!(y_array,x)     #testing Base.copyto!(dest::Array,A::ArrayPartition)