/
projectors.jl
110 lines (93 loc) · 2.9 KB
/
projectors.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
"""
Projector(ϕs, orbitals, S)
Represents the projector on the subspace spanned by the radial
orbitals `ϕs` (corresponding to `orbitals`).
"""
struct Projector{T,B<:AbstractQuasiMatrix,RO<:RadialOrbital{T,B},O,Metric,Selection} <: NBodyOperator{1}
ϕs::Vector{RO}
orbitals::Vector{O}
S::Metric
sel::Selection
function Projector(ϕs::Vector{RO}, orbitals::Vector{O}, S::Metric, sel=nothing) where {T,B,RO<:RadialOrbital{T,B},O,Metric}
m,n = size(S)
all(ϕ -> length(ϕ.args[2]) == m, ϕs) ||
throw(DimensionMismatch())
if isnothing(sel)
sel = m == n ? Colon() : 1:m
end
new{T,B,RO,O,Metric,typeof(sel)}(ϕs, orbitals, S, sel)
end
end
Base.eltype(projector::Projector{T}) where T = T
function Base.size(projector::Projector)
n = size(projector.S,2)
n,n
end
Base.size(projector::Projector, i) = size(projector)[i]
Base.similar(projector::Projector{T}) where T = Matrix{T}(undef, size(projector))
Base.iszero(me::OrbitalMatrixElement{1,A,<:Projector,B}) where {A<:SpinOrbital,B<:SpinOrbital} =
me.a[1] ∉ me.o.orbitals || me.b[1] ∉ me.o.orbitals
function Base.show(io::IO, projector::Projector)
write(io, "P(")
write(io, join(string.(projector.orbitals), " "))
write(io, ")")
if !(projector.sel isa Colon)
write(io, " restricted to elements $(projector.sel)")
end
end
# * Multiplication by projector
proj_view(y, ::Colon) = y
proj_view(y, sel) = view(y, sel)
function LinearAlgebra.mul!(y::AbstractVector, p::Projector, x::AbstractVector,
α::Number=true, β::Number=false)
if iszero(β)
y .= false
elseif !isone(β)
lmul!(β, y)
end
iszero(α) && return y
yv = proj_view(y, p.sel)
for ϕ in p.ϕs
ϕv = ϕ.args[2]
c = α*dot(ϕv, p.S, x)
yv .+= c .* ϕv
end
y
end
function LinearAlgebra.mul!(Y::AbstractMatrix, p::Projector, X::AbstractMatrix,
α::Number=true, β::Number=false)
n = size(Y,2)
for j = 1:n
mul!(view(Y, :, j), p, view(X, :, j), α, β)
end
Y
end
LinearAlgebra.Matrix(p::Projector) =
mul!(similar(p), p, I(size(p.S,2)))
# * Rejection
projectout!(y::RO, ::Nothing) where RO = y
function projectout!(y::AbstractVector{T}, projector::Projector) where T
s = zero(T)
yv = proj_view(y, projector.sel)
for ϕ in projector.ϕs
c = dot(ϕ.args[2], projector.S, y)
yv .-= c .* ϕ.args[2]
s += c
end
s
end
function projectout!(Y::AbstractMatrix{T}, projector::Projector) where T
s = zero(T)
n = size(Y, 2)
for j = 1:n
s += projectout!(view(Y, :, j), projector)
end
s
end
"""
projectout!(y, projector)
Project out all components of `y` parallel to the radial orbitals
`projector.ϕs`.
"""
projectout!(y::RadialOrbital, projector::Projector) =
projectout!(y.args[2], projector)