Make alloc_matmul type-stable

src/FilterMD.jl

@@ -32,10 +32,11 @@ dimension ordering. In many cases, this algorithm exhibits the best cache behavi
 """ ->
 A_mul_B_perm(M::AbstractMatrix, src, dim::Integer) = A_mul_B_perm!(alloc_matmul(M,src,dim), M, src, dim)
 
-function alloc_matmul(M,src,dim)
+function alloc_matmul{S,N}(M,src::AbstractArray{S,N},dim)
     sz = [size(src)...]
     sz[dim] = size(M,1)
-    Array(promote_type(eltype(M), eltype(src)), sz...)
+    T = promote_type(eltype(M), S)
+    Array(T, sz...)::Array{T,N}
 end
 
 include("tridiag.jl")

src/woodbury.jl

@@ -5,10 +5,14 @@ function _A_ldiv_B_md!(dest, W::Woodbury, src,  R1, R2)
     tmp3 = _A_mul_B_md(W.U, tmp2, R1, R2)
     # TODO?: would be nice to fuse the next two steps
     tmp4 = _A_ldiv_B_md(W.A, tmp3, R1, R2)
-    for I in eachindex(dest, tmp4)
-        dest[I] -= tmp4[I]
+    sub!(dest, tmp4)
+end
+
+function sub!(A, B)
+    for I in eachindex(A, B)
+        A[I] -= B[I]
     end
-    dest
+    A
 end
 
 check_matrix(W::Woodbury) = check_matrix(W.A)