From e6e72aedaa0d779e27bb525314fe428c25d35060 Mon Sep 17 00:00:00 2001 From: mtanneau Date: Thu, 22 Oct 2020 20:33:14 -0400 Subject: [PATCH] Pardiso linear solver --- Project.toml | 1 + src/KKT/KKT.jl | 1 + src/KKT/pardiso.jl | 183 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 185 insertions(+) create mode 100644 src/KKT/pardiso.jl diff --git a/Project.toml b/Project.toml index c84bd4af..3c6cd737 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QPSReader = "10f199a5-22af-520b-b891-7ce84a7b1bd0" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl index f193cc2a..aa1c1433 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -108,6 +108,7 @@ include("lapack.jl") include("cholmod.jl") include("ldlfact.jl") include("krylov.jl") +include("pardiso.jl") """ default_options(T) diff --git a/src/KKT/pardiso.jl b/src/KKT/pardiso.jl new file mode 100644 index 00000000..399a8ca4 --- /dev/null +++ b/src/KKT/pardiso.jl @@ -0,0 +1,183 @@ +import Pardiso + +mutable struct MKLPardisoSQD <: AbstractKKTSolver{Float64} + m::Int # Number of rows + n::Int # Number of columns + + # Problem data + A::SparseMatrixCSC{Float64, Int} + θ::Vector{Float64} + regP::Vector{Float64} # primal regularization + regD::Vector{Float64} # dual regularization + + # Left-hand side matrix + S::SparseMatrixCSC{Float64, Int} + + # Linear solver + ps::Pardiso.MKLPardisoSolver + + function MKLPardisoSQD(A::SparseMatrixCSC{Float64}) + + m, n = size(A) + θ = ones(n) + + # We store we lower-triangular of the matrix + S = [ + spdiagm(0 => -θ) spzeros(n, m); + A spdiagm(0 => ones(m)) + ] + + ps = Pardiso.MKLPardisoSolver() + + # We use symmetric indefinite matrices + Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_INDEF) + Pardiso.pardisoinit(ps) + + # Set number of threads + Pardiso.set_nprocs!(ps, 1) + + # Do the analysis + Pardiso.set_phase!(ps, Pardiso.ANALYSIS) + Pardiso.pardiso(ps, S, ones(m+n)) + + return new(m, n, A, θ, ones(Float64, n), ones(Float64, m), S, ps) + + return kkt + end +end + +mutable struct PardisoSQD <: AbstractKKTSolver{Float64} + m::Int # Number of rows + n::Int # Number of columns + + # Problem data + A::SparseMatrixCSC{Float64, Int} + θ::Vector{Float64} + regP::Vector{Float64} # primal regularization + regD::Vector{Float64} # dual regularization + + # Left-hand side matrix + S::SparseMatrixCSC{Float64, Int} + + # Linear solver + ps::Pardiso.PardisoSolver + + function PardisoSQD(A::SparseMatrixCSC{Float64}) + + m, n = size(A) + θ = ones(n) + + # We store we lower-triangular of the matrix + S = [ + spdiagm(0 => -θ) A'; + A spdiagm(0 => ones(m)) + ] + + ps = Pardiso.PardisoSolver() + + # We use symmetric indefinite matrices + Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_INDEF) + Pardiso.pardisoinit(ps) + + S_pardiso = Pardiso.get_matrix(ps, S, :N) + @assert istril(S_pardiso) + + # Use a direct method + Pardiso.set_solver!(ps, 1) + + # Do the analysis + Pardiso.set_phase!(ps, Pardiso.ANALYSIS) + Pardiso.pardiso(ps, S_pardiso, ones(m+n)) + + return new(m, n, A, θ, ones(Float64, n), ones(Float64, m), S_pardiso, ps) + + return kkt + end +end + +setup(::Type{MKLPardisoSQD}, A) = MKLPardisoSQD(A) +backend(::MKLPardisoSQD) = "MKLPardiso" +linear_system(::MKLPardisoSQD) = "Augmented system" + +setup(::Type{PardisoSQD}, A) = PardisoSQD(A) +backend(::PardisoSQD) = "Pardiso" +linear_system(::PardisoSQD) = "Augmented system" + +""" + update!(kkt, θ, regP, regD) + +Update LDLᵀ factorization of the augmented system. + +Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute + the factorization. +Throws a `PosDefException` if matrix is not quasi-definite. +""" +function update!( + kkt::Union{MKLPardisoSQD, PardisoSQD}, + θ::AbstractVector{Float64}, + regP::AbstractVector{Float64}, + regD::AbstractVector{Float64} +) + # Sanity checks + length(θ) == kkt.n || throw(DimensionMismatch( + "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." + )) + length(regP) == kkt.n || throw(DimensionMismatch( + "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" + )) + length(regD) == kkt.m || throw(DimensionMismatch( + "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" + )) + + m, n = kkt.m, kkt.n + + # Update diagonal scaling + kkt.θ .= θ + # Update regularizers + kkt.regP .= regP + kkt.regD .= regD + + # Update S. + # S is stored as lower-triangular and only its diagonal changes. + @inbounds for j in 1:kkt.n + k = kkt.S.colptr[j] + kkt.S.nzval[k] = -kkt.θ[j] - regP[j] + end + @inbounds for i in 1:kkt.m + k = kkt.S.colptr[kkt.n+i] + kkt.S.nzval[k] = regD[i] + end + + # Compute numerical factorization + Pardiso.set_phase!(kkt.ps, Pardiso.NUM_FACT) + Pardiso.pardiso(kkt.ps, kkt.S, zeros(kkt.m + kkt.n)) + + return nothing +end + +""" + solve!(dx, dy, kkt, ξp, ξd) + +Solve the augmented system, overwriting `dx, dy` with the result. +""" +function solve!( + dx::Vector{Float64}, dy::Vector{Float64}, + kkt::Union{MKLPardisoSQD, PardisoSQD}, + ξp::Vector{Float64}, ξd::Vector{Float64} +) + m, n = kkt.m, kkt.n + + # Set-up right-hand side + ξ = [ξd; ξp] + + # Solve augmented system + d = zeros(m + n) + Pardiso.set_phase!(kkt.ps, Pardiso.SOLVE_ITERATIVE_REFINE) + Pardiso.pardiso(kkt.ps, d, kkt.S, ξ) + + # Recover dx, dy + @views dx .= d[1:n] + @views dy .= d[(n+1):(m+n)] + + return nothing +end \ No newline at end of file