Skip to content

Commit 490284e

Browse files
dlfivefiftydlfivefifty
committed
Add reflector!, Apply
1 parent 669bf08 commit 490284e

File tree

4 files changed

+184
-3
lines changed

4 files changed

+184
-3
lines changed

‎src/ArrayLayouts.jl

Lines changed: 65 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -46,11 +46,11 @@ else
4646
import Base: require_one_based_indexing
4747
end
4848

49-
export materialize, materialize!, MulAdd, muladd!, Ldiv, Lmul, Rmul, MemoryLayout, AbstractStridedLayout,
49+
export materialize, materialize!, MulAdd, muladd!, Ldiv, Lmul, Rmul, lmul, rmul, MemoryLayout, AbstractStridedLayout,
5050
DenseColumnMajor, ColumnMajor, ZerosLayout, FillLayout, AbstractColumnMajor, RowMajor, AbstractRowMajor,
5151
DiagonalLayout, ScalarLayout, SymTridiagonalLayout, HermitianLayout, SymmetricLayout, TriangularLayout,
5252
UnknownLayout, AbstractBandedLayout, ApplyBroadcastStyle, ConjLayout, AbstractFillLayout,
53-
colsupport, rowsupport, lazy_getindex
53+
colsupport, rowsupport, lazy_getindex, QLayout
5454

5555
struct ApplyBroadcastStyle <: BroadcastStyle end
5656
@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle})
@@ -64,12 +64,75 @@ include("lmul.jl")
6464
include("ldiv.jl")
6565
include("diagonal.jl")
6666
include("triangular.jl")
67+
include("factorizations.jl")
6768

6869
@inline sub_materialize(_, V) = Array(V)
6970
@inline sub_materialize(V::SubArray) = sub_materialize(MemoryLayout(typeof(V)), V)
7071

7172
@inline lazy_getindex(A, I...) = sub_materialize(view(A, I...))
7273

74+
zero!(A::AbstractArray{T}) where T = fill!(A,zero(T))
75+
function zero!(A::AbstractArray{<:AbstractArray})
76+
for a in A
77+
zero!(a)
78+
end
79+
A
80+
end
81+
7382
_fill_lmul!(β, A::AbstractArray{T}) where T = iszero(β) ? zero!(A) : lmul!(β, A)
7483

84+
85+
# Elementary reflection similar to LAPACK. The reflector is not Hermitian but
86+
# ensures that tridiagonalization of Hermitian matrices become real. See lawn72
87+
@inline function reflector!(x::AbstractVector)
88+
require_one_based_indexing(x)
89+
n = length(x)
90+
n == 0 && return zero(eltype(x))
91+
@inbounds begin
92+
ξ1 = x[1]
93+
normu = abs2(ξ1)
94+
for i = 2:n
95+
normu += abs2(x[i])
96+
end
97+
if iszero(normu)
98+
return zero(ξ1/normu)
99+
end
100+
normu = sqrt(normu)
101+
ν = copysign(normu, real(ξ1))
102+
ξ1 += ν
103+
x[1] = -ν
104+
for i = 2:n
105+
x[i] /= ξ1
106+
end
107+
end
108+
ξ1/ν
109+
end
110+
111+
# apply reflector from left
112+
@inline function reflectorApply!(x::AbstractVector, τ::Number, A::AbstractVecOrMat)
113+
m,n = size(A,1),size(A,2)
114+
if length(x) != m
115+
throw(DimensionMismatch("reflector has length $(length(x)), which must match the first dimension of matrix A, $m"))
116+
end
117+
m == 0 && return A
118+
@inbounds begin
119+
for j = 1:n
120+
# dot
121+
vAj = A[1, j]
122+
for i = 2:m
123+
vAj += x[i]'*A[i, j]
124+
end
125+
126+
vAj = conj(τ)*vAj
127+
128+
# ger
129+
A[1, j] -= vAj
130+
for i = 2:m
131+
A[i, j] -= x[i]*vAj
132+
end
133+
end
134+
end
135+
return A
136+
end
137+
75138
end

‎src/factorizations.jl

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
struct QLayout <: MemoryLayout end
2+
3+
MemoryLayout(::Type{<:AbstractQ}) = QLayout()
4+
5+
transposelayout(::QLayout) = QLayout()
6+
7+
8+
copy(M::Lmul{QLayout}) = copyto!(similar(M), M)
9+
10+
function copyto!(dest::AbstractArray{T}, M::Lmul{QLayout}) where T
11+
A,B = M.A,M.B
12+
if size(dest,1) == size(B,1)
13+
copyto!(dest, B)
14+
else
15+
copyto!(view(dest,1:size(B,1),:), B)
16+
zero!(@view(dest[size(B,1)+1:end,:]))
17+
end
18+
materialize!(Lmul(A,dest))
19+
end
20+
21+
function copyto!(dest::AbstractArray, M::Ldiv{QLayout})
22+
A,B = M.A,M.B
23+
copyto!(dest, B)
24+
materialize!(Ldiv(A,dest))
25+
end
26+
27+
materialize!(M::Ldiv{QLayout}) = materialize!(Lmul(M.A',M.B))

‎src/lmul.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,9 @@ const BlasMatRmulMat{StyleA,StyleB,T<:BlasFloat} = Rmul{StyleA,StyleB,<:Abstract
4141

4242
axes(M::MatLmulVec) = (axes(M.A,1),)
4343

44+
lmul(A, B) = materialize(Lmul(A, B))
45+
46+
materialize(L::Lmul) = copy(instantiate(L))
4447

4548
copy(M::Lmul) = materialize!(Lmul(M.A,copy(M.B)))
4649
copy(M::Rmul) = materialize!(Rmul(copy(M.A),M.B))

‎src/muladd.jl

Lines changed: 89 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,13 @@ function check_mul_axes(A, B, C...)
4949
check_mul_axes(B, C...)
5050
end
5151

52+
# we need to special case AbstractQ as it allows non-compatiple multiplication
53+
function check_mul_axes(A::AbstractQ, B, C...)
54+
axes(A.factors, 1) == axes(B, 1) || axes(A.factors, 2) == axes(B, 1) ||
55+
throw(DimensionMismatch("First axis of B, $(axes(B,1)) must match either axes of A, $(axes(A))"))
56+
check_mul_axes(B, C...)
57+
end
58+
5259

5360
function instantiate(M::MulAdd)
5461
@boundscheck check_mul_axes(M.α, M.A, M.B)
@@ -373,4 +380,85 @@ copy(M::MulAdd{<:AbstractFillLayout,<:AbstractFillLayout,<:AbstractFillLayout})
373380
copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.B.diag)) .* M.A .+ M.β .* M.C
374381
copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout},ZerosLayout}) = (M.α * getindex_value(M.B.diag)) .* M.A
375382

376-
BroadcastStyle(::Type{<:MulAdd}) = ApplyBroadcastStyle()
383+
BroadcastStyle(::Type{<:MulAdd}) = ApplyBroadcastStyle()
384+
385+
scalarone(::Type{T}) where T = one(T)
386+
scalarone(::Type{<:AbstractArray{T}}) where T = scalarone(T)
387+
scalarzero(::Type{T}) where T = zero(T)
388+
scalarzero(::Type{<:AbstractArray{T}}) where T = scalarzero(T)
389+
390+
fillzeros(::Type{T}, ax) where T = Zeros{T}(ax)
391+
392+
function mul!(dest::AbstractArray{W}, A::AbstractArray{T}, b::AbstractArray{V}) where {T,V,W}
393+
TVW = promote_type(W, _mul_eltype(T,V))
394+
muladd!(scalarone(TVW), A, b, scalarzero(TVW), dest)
395+
end
396+
397+
function MulAdd(A::AbstractArray{T}, B::AbstractVector{V}) where {T,V}
398+
TV = _mul_eltype(eltype(A), eltype(B))
399+
MulAdd(scalarone(TV), A, B, scalarzero(TV), fillzeros(TV,(axes(A,1))))
400+
end
401+
402+
function MulAdd(A::AbstractArray{T}, B::AbstractMatrix{V}) where {T,V}
403+
TV = _mul_eltype(eltype(A), eltype(B))
404+
MulAdd(scalarone(TV), A, B, scalarzero(TV), fillzeros(TV,(axes(A,1),axes(B,2))))
405+
end
406+
407+
mul(A::AbstractArray, B::AbstractArray) = materialize(MulAdd(A,B))
408+
409+
macro lazymul(Typ)
410+
ret = quote
411+
LinearAlgebra.mul!(dest::AbstractVector, A::$Typ, b::AbstractVector) =
412+
ArrayLayouts.mul!(dest,A,b)
413+
414+
LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::AbstractMatrix) =
415+
ArrayLayouts.mul!(dest,A,b)
416+
LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Typ) =
417+
ArrayLayouts.mul!(dest,A,b)
418+
419+
Base.:*(A::$Typ, B::$Typ) = ArrayLayouts.mul(A,B)
420+
Base.:*(A::$Typ, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
421+
Base.:*(A::$Typ, B::AbstractVector) = ArrayLayouts.mul(A,B)
422+
Base.:*(A::AbstractMatrix, B::$Typ) = ArrayLayouts.mul(A,B)
423+
Base.:*(A::LinearAlgebra.AdjointAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
424+
Base.:*(A::LinearAlgebra.TransposeAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
425+
426+
Base.:*(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.lmul(A,B)
427+
Base.:*(A::$Typ, B::LinearAlgebra.AbstractQ) = ArrayLayouts.rmul(A,B)
428+
end
429+
for Struc in (:AbstractTriangular, :Diagonal)
430+
ret = quote
431+
$ret
432+
433+
Base.:*(A::LinearAlgebra.$Struc, B::$Typ) = ArrayLayouts.mul(A,B)
434+
Base.:*(A::$Typ, B::LinearAlgebra.$Struc) = ArrayLayouts.mul(A,B)
435+
end
436+
end
437+
for Mod in (:Adjoint, :Transpose, :Symmetric, :Hermitian)
438+
ret = quote
439+
$ret
440+
441+
LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Mod{<:Any,<:AbstractMatrix}) =
442+
ArrayLayouts.mul!(dest,A,b)
443+
444+
LinearAlgebra.mul!(dest::AbstractVector, A::$Mod{<:Any,<:$Typ}, b::AbstractVector) =
445+
ArrayLayouts.mul!(dest,A,b)
446+
447+
Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
448+
Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
449+
Base.:*(A::AbstractMatrix, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
450+
Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractVector) = ArrayLayouts.mul(A,B)
451+
452+
Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.mul(A,B)
453+
Base.:*(A::$Typ, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
454+
455+
Base.:*(A::$Mod{<:Any,<:$Typ}, B::Diagonal) = ArrayLayouts.mul(A,B)
456+
Base.:*(A::Diagonal, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
457+
458+
Base.:*(A::LinearAlgebra.AbstractTriangular, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
459+
Base.:*(A::$Mod{<:Any,<:$Typ}, B::LinearAlgebra.AbstractTriangular) = ArrayLayouts.mul(A,B)
460+
end
461+
end
462+
463+
esc(ret)
464+
end

0 commit comments

Comments
 (0)