Static arrays (#99)

* Initial changes for StaticArrays [ci skip]

* Changes to reduce allocation

* Refactoring to reduce allocations

* Slight changes in results

* Relax tolerance on comparison

Author: dmbates

REQUIRE

@@ -8,5 +8,6 @@ GLM 0.7
 NamedArrays
 NLopt 0.3
 Showoff
+StaticArrays 0.6
 StatsBase 0.11
 StatsFuns 0.3

src/MixedModels.jl

@@ -3,7 +3,7 @@ __precompile__()
 module MixedModels
 
 using ArgCheck, BlockArrays, DataArrays, DataFrames
-using Distributions, GLM, NLopt, Showoff, StatsBase
+using Distributions, GLM, NLopt, Showoff, StaticArrays, StatsBase
 using StatsFuns: log2π
 using NamedArrays: NamedArray, setnames!
 using Base.LinAlg: BlasFloat, BlasReal, HermOrSym, PosDefException, checksquare, copytri!

src/linalg.jl

@@ -25,45 +25,12 @@ function αβA_mul_Bc!(α::T, A::SparseMatrixCSC{T}, B::SparseMatrixCSC{T},
     C
 end
 
-#= This is not tested.  See if it is really needed.
-function αβA_mul_Bc!(α::T, A::SparseMatrixCSC{T}, B::SparseMatrixCSC{T},
-                     β::T, C::SparseMatrixCSC{T}) where T <: Number
-    @argcheck B.m == C.n && A.m == C.m && A.n == B.n  DimensionMismatch
-    anz = nonzeros(A)
-    arv = rowvals(A)
-    bnz = nonzeros(B)
-    brv = rowvals(B)
-    cnz = nonzeros(C)
-    crv = rowvals(C)
-    if β ≠ one(T)
-        iszero(β) ? fill!(cnz, β) : scale!(cnz, β)
-    end
-    for j = 1:A.n
-        for ib in nzrange(B, j)
-            αbnz = α * bnz[ib]
-            jj = brv[ib]
-            for ia in nzrange(A, j)
-                crng = nzrange(C, jj)
-                ind = findfirst(crv[crng], arv[ia])
-                if iszero(ind)
-                    throw(ArgumentError("A*B' has nonzero positions not in C"))
-                end
-                cnz[crng[ind]] += anz[ia] * αbnz
-            end
-        end
-    end
-    C
-end
-=#
-
 function αβA_mul_Bc!(α::T, A::StridedVecOrMat{T}, B::SparseMatrixCSC{T}, β::T,
                      C::StridedVecOrMat{T}) where T
     m, n = size(A)
     p, q = size(B)
     r, s = size(C)
-    if r ≠ m || s ≠ p || n ≠ q
-        throw(DimensionMismatch("size(C,1) ≠ size(A,1) or size(C,2) ≠ size(B,1) or size(A,2) ≠ size(B,2)"))
-    end
+    @argcheck(r == m && s == p && n == q, DimensionMismatch)
     if β ≠ one(T)
         iszero(β) ? fill!(C, β) : scale!(C, β)
     end

src/linalg/cholUnblocked.jl

@@ -29,4 +29,9 @@ function cholUnblocked!(A::StridedMatrix{T}, ::Type{Val{:L}}) where T<:BlasFloat
     A
 end
 
-cholUnblocked!(D::UniformBlockDiagonal, ::Type{Val{:L}}) = cholUnblocked!.(D.facevec, Val{:L})
+function cholUnblocked!(D::UniformBlockDiagonal, ::Type{Val{:L}})
+    for f in D.facevec
+        cholUnblocked!(f, Val{:L})
+    end
+    D
+end

src/linalg/logdet.jl

@@ -32,10 +32,11 @@ end
 Return the value of `log(det(Λ'Z'ZΛ + I))` evaluated in place.
 """
 function logdet(m::LinearMixedModel{T}) where {T}
-    s = zero(T)
+    s = log(one(T))
+    Ldat = m.L.data
     for (i, trm) in enumerate(m.trms)
         if isa(trm, AbstractFactorReTerm)
-            s += T(LD(m.L.data[Block(i, i)]))
+            s += LD(m.L.data[Block(i, i)])
         end
     end
     2s

src/linalg/rankUpdate.jl

@@ -15,17 +15,22 @@ function rankUpdate!(α::T, a::StridedVector{T},
     BLAS.syr!(A.uplo, α, a, A.data)
     A
 end
-rankUpdate!(a::StridedVector{T}, A::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix} = rankUpdate!(one(T), a, A)
+
+rankUpdate!(a::StridedVector{T}, A::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix} =
+    rankUpdate!(one(T), a, A)
 
 rankUpdate!(α::T, A::StridedMatrix{T}, β::T,
-C::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix} = BLAS.syrk!(C.uplo, 'N', α, A, β, C.data)
+            C::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix} =
+    BLAS.syrk!(C.uplo, 'N', α, A, β, C.data)
+
 rankUpdate!(α::T, A::StridedMatrix{T}, C::HermOrSym{T,S}) where {T<:Real,S<:StridedMatrix} =
     rankUpdate!(α, A, one(T), C)
+
 rankUpdate!(A::StridedMatrix{T}, C::HermOrSym{T,S}) where {T<:Real,S<:StridedMatrix} =
     rankUpdate!(one(T), A, one(T), C)
 
-function rankUpdate!(α::T, A::SparseMatrixCSC{T,I},
-                     β::T, C::HermOrSym{T,S}) where {T,I,S<:StridedMatrix{T}}
+function rankUpdate!(α::T, A::SparseMatrixCSC{T},
+                     β::T, C::HermOrSym{T,S}) where {T,S<:StridedMatrix{T}}
     m, n = size(A)
     @argcheck m == size(C, 2) && C.uplo == 'L' DimensionMismatch
     Cd = C.data
@@ -73,6 +78,7 @@ function rankUpdate!(α::T, A::SparseMatrixCSC{T},
                      C::HermOrSym{T,UniformBlockDiagonal{T}}) where T<:Number
     m, n, k = size(C.data.data)
     @argcheck m == n && size(A, 1) == m * k DimensionMismatch
+    # Another expensive evaluation in terms of storage allocation
     aat = α * (A * A')
     nz = nonzeros(aat)
     rv = rowvals(aat)