Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

WIP: start using generic 5-arg mul! #308

Merged
merged 4 commits into from
Aug 15, 2019
Merged

WIP: start using generic 5-arg mul! #308

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
74630e9
start using generic 5-arg mul! in some places
chriscoey Aug 15, 2019
11e31b8
use 5-arg mul! for HypBlockMatrix
chriscoey Aug 15, 2019
61866d8
start simplifying ops in qrchol
chriscoey Aug 15, 2019
8ac8005
further simplify qrchol
chriscoey Aug 15, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 12 additions & 16 deletions src/Solvers/homogeneous_self_dual/qrchol.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,10 +41,9 @@ mutable struct QRCholCombinedHSDSystemSolver{T <: HypReal} <: CombinedHSDSystemS
GQ2
QpbxGHbz
Q1pbxGHbz
Q2pbxGHbz
Q2div
GQ1x
HGQ1x
Q2div
HGQ2
Q2GHGQ2
Gxi
Expand Down Expand Up @@ -118,10 +117,9 @@ mutable struct QRCholCombinedHSDSystemSolver{T <: HypReal} <: CombinedHSDSystemS
end
system_solver.QpbxGHbz = Matrix{T}(undef, n, 3)
system_solver.Q1pbxGHbz = view(system_solver.QpbxGHbz, 1:p, :)
system_solver.Q2pbxGHbz = view(system_solver.QpbxGHbz, (p + 1):n, :)
system_solver.Q2div = view(system_solver.QpbxGHbz, (p + 1):n, :)
system_solver.GQ1x = Matrix{T}(undef, q, 3)
system_solver.HGQ1x = similar(system_solver.GQ1x)
system_solver.Q2div = Matrix{T}(undef, nmp, 3)
system_solver.Gxi = similar(system_solver.GQ1x)
system_solver.HGxi = similar(system_solver.Gxi)

Expand Down Expand Up @@ -173,10 +171,9 @@ function get_combined_directions(solver::HSDSolver{T}, system_solver::QRCholComb
GQ2 = system_solver.GQ2
QpbxGHbz = system_solver.QpbxGHbz
Q1pbxGHbz = system_solver.Q1pbxGHbz
Q2pbxGHbz = system_solver.Q2pbxGHbz
Q2div = system_solver.Q2div
GQ1x = system_solver.GQ1x
HGQ1x = system_solver.HGQ1x
Q2div = system_solver.Q2div
HGQ2 = system_solver.HGQ2
Q2GHGQ2 = system_solver.Q2GHGQ2
Gxi = system_solver.Gxi
Expand Down Expand Up @@ -233,15 +230,14 @@ function get_combined_directions(solver::HSDSolver{T}, system_solver::QRCholComb

ldiv!(model.Ap_R', yi)

mul!(QpbxGHbz, model.G', zi)
@. QpbxGHbz += xi
copyto!(QpbxGHbz, xi)
mul!(QpbxGHbz, model.G', zi, true, true)
lmul!(model.Ap_Q', QpbxGHbz)

if !iszero(size(Q2div, 1))
mul!(GQ1x, GQ1, yi)
block_hessian_product!(HGQ1x_k, GQ1x_k)
mul!(Q2div, GQ2', HGQ1x)
@. Q2div = Q2pbxGHbz - Q2div
mul!(Q2div, GQ2', HGQ1x, -1, true)

block_hessian_product!(HGQ2_k, GQ2_k)
mul!(Q2GHGQ2, GQ2', HGQ2)
Expand Down Expand Up @@ -280,8 +276,8 @@ function get_combined_directions(solver::HSDSolver{T}, system_solver::QRCholComb
end
end

@. xi1 = yi
@. xi2 = Q2div
copyto!(xi1, yi)
copyto!(xi2, Q2div)
lmul!(model.Ap_Q, xi)

mul!(Gxi, model.G, xi)
Expand All @@ -290,8 +286,8 @@ function get_combined_directions(solver::HSDSolver{T}, system_solver::QRCholComb
@. zi = HGxi - zi

if !iszero(length(yi))
mul!(yi, GQ1', HGxi)
@. yi = Q1pbxGHbz - yi
copyto!(yi, Q1pbxGHbz)
mul!(yi, GQ1', HGxi, -1, true)
ldiv!(model.Ap_R, yi)
end

Expand All @@ -303,8 +299,8 @@ function get_combined_directions(solver::HSDSolver{T}, system_solver::QRCholComb
@. x += tau_sol * x1
@. y += tau_sol * y1
@. z += tau_sol * z1
mul!(s, model.G, x)
@. s = -s + tau_sol * model.h
copyto!(s, model.h)
mul!(s, model.G, x, -1, tau_sol)
kap_sol = -dot(model.c, x) - dot(model.b, y) - dot(model.h, z) - tau_rhs
return (tau_sol, kap_sol)
end
Expand Down
5 changes: 3 additions & 2 deletions src/Solvers/homogeneous_self_dual/solver.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,10 +198,11 @@ function calc_residual(solver::HSDSolver{T}) where {T <: HypReal}
# x_residual = -A'*y - G'*z - c*tau
x_residual = solver.x_residual
mul!(x_residual, model.G', point.z)
x_residual .= -model.A' * point.y - x_residual # TODO remove allocs
mul!(x_residual, model.A', point.y, true, true)
solver.x_norm_res_t = norm(x_residual)
@. x_residual -= model.c * solver.tau
@. x_residual += model.c * solver.tau
solver.x_norm_res = norm(x_residual) / solver.tau
@. x_residual *= -1

# y_residual = A*x - b*tau
y_residual = solver.y_residual
Expand Down
46 changes: 5 additions & 41 deletions src/linearalgebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,61 +52,25 @@ size(A::HypBlockMatrix, d) = (d == 1 ? A.nrows : A.ncols)

adjoint(A::HypBlockMatrix{T}) where {T <: HypReal} = HypBlockMatrix{T}(A.ncols, A.nrows, adjoint.(A.blocks), A.cols, A.rows)

function mul!(y::AbstractVector{T}, A::HypBlockMatrix{T}, x::AbstractVector{T}) where {T <: HypReal}
# TODO try to speed up by using better logic for alpha and beta (see Julia's 5-arg mul! code)
# TODO check that this eliminates allocs when using IterativeSolvers methods, and that it is as fast as possible
function mul!(y::AbstractVector{T}, A::HypBlockMatrix{T}, x::AbstractVector{T}, alpha::Number, beta::Number) where {T <: HypReal}
@assert size(x, 1) == A.ncols
@assert size(y, 1) == A.nrows
@assert size(x, 2) == size(y, 2)

y .= zero(T)

@. y *= beta
for (b, r, c) in zip(A.blocks, A.rows, A.cols)
if isempty(r) || isempty(c)
continue
end
# println()
# if b isa UniformScaling
# println("I")
# else
# println(size(b))
# end
# println(r, " , ", c)
xk = view(x, c)
yk = view(y, r)
yk .+= b * xk # TODO need inplace mul+add
# mul!(yk, b, xk, α = one(T), β = one(T))
mul!(yk, b, xk, alpha, true)
end

return y
end

# function mul!(y::AbstractVector{T}, A::Adjoint{T, HypBlockMatrix{T}}, x::AbstractVector{T}) where {T <: HypReal}
# @assert size(x, 1) == A.ncols
# @assert size(y, 1) == A.nrows
# @assert size(x, 2) == size(y, 2)
#
# y .= zero(T)
#
# for (b, r, c) in zip(A.blocks, A.rows, A.cols)
# if isempty(r) || isempty(c)
# continue
# end
# # println()
# # if b isa UniformScaling
# # println("I")
# # else
# # println(size(b))
# # end
# # println(r, " , ", c)
# xk = view(x, r)
# yk = view(y, c)
# yk .+= b' * xk # TODO need inplace mul+add
# # mul!(yk, b', xk)
# # mul!(yk, b', xk, α = one(T), β = one(T))
# end
#
# return y
# end

*(A::HypBlockMatrix{T}, x::AbstractVector{T}) where {T <: HypReal} = mul!(similar(x, size(A, 1)), A, x)

-(A::HypBlockMatrix{T}) where {T <: HypReal} = HypBlockMatrix{T}(A.nrows, A.ncols, -A.blocks, A.rows, A.cols)