Skip to content

Commit a560d35

Browse files
committed
Changes to reduce allocation
1 parent db71918 commit a560d35

File tree

3 files changed

+26
-12
lines changed

3 files changed

+26
-12
lines changed

src/linalg.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -126,9 +126,9 @@ function A_rdiv_Bc!(A::Matrix{T}, B::LowerTriangular{T,UniformBlockDiagonal{T}})
126126
m, n, k = size(B.data.data)
127127
@argcheck size(A, 2) == size(B, 1) && m == n DimensionMismatch
128128
offset = 0
129+
one2m = 1:m
129130
for f in B.data.facevec
130-
## FIXME call BLAS.trsm directly
131-
A_rdiv_Bc!(view(A, :, (1:m) + offset), LowerTriangular(f))
131+
BLAS.trsm!('R', 'L', 'T', 'N', one(T), f, view(A, :, one2m + offset))
132132
offset += m
133133
end
134134
A

src/linalg/logdet.jl

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -8,12 +8,23 @@ Return `log(det(tril(A)))` evaluated in place.
88
LD(d::Diagonal{T}) where {T<:Number} = sum(log, d.diag)
99

1010
function LD(d::UniformBlockDiagonal{T}) where T
11-
m, n, k = size(d.data)
12-
dind = diagind(m, n)
13-
sum(log, f[i] for f in d.facevec, i in dind)
11+
dat = d.data
12+
m, n, k = size(dat)
13+
m == n || throw(ArgumentError("Blocks of d must be square"))
14+
s = log(one(T))
15+
@inbounds for j in 1:k, i in 1:m
16+
s += log(dat[i,i,j])
17+
end
18+
s
1419
end
1520

16-
LD(d::DenseMatrix{T}) where {T} = sum(log, d[i] for i in diagind(d))
21+
function LD(d::DenseMatrix{T}) where T
22+
s = log(one(T))
23+
for i in 1:Base.LinAlg.checksquare(d)
24+
s += log(d[i, i])
25+
end
26+
s
27+
end
1728

1829
"""
1930
logdet(m::LinearMixedModel)

src/linalg/scaleInflate.jl

Lines changed: 9 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -36,13 +36,16 @@ function scaleInflate!(Ljj::LowerTriangular{T,UniformBlockDiagonal{T}},
3636
Ajj::UniformBlockDiagonal{T},
3737
Λj::VectorFactorReTerm{T}) where {T}
3838
@argcheck size(Ljj) == size(Ajj) DimensionMismatch
39+
Ljjdat = Ljj.data
40+
Ljjdd = Ljjdat.data
41+
copy!(Ljjdd, Ajj.data)
42+
k, m, n = size(Ljjdd)
3943
λ = LowerTriangular(Λj.Λ)
40-
k = vsize(Λj)
41-
for (Lf, Af) in zip(Ljj.data.facevec, Ajj.facevec)
42-
Ac_mul_B!(λ, A_mul_B!(copy!(Lf, Af), λ))
43-
for j in 1:k
44-
Lf[j, j] += one(T)
45-
end
44+
for Lf in Ljjdat.facevec
45+
Ac_mul_B!(λ, A_mul_B!(Lf, λ))
46+
end
47+
for j in 1:n, i in 1:k
48+
Ljjdd[i, i, j] += one(T)
4649
end
4750
Ljj
4851
end

0 commit comments

Comments
 (0)