diff --git a/Project.toml b/Project.toml index 768294c6..4bf4c6bc 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1" KLU = "ef3ab10e-7fda-4108-b977-705223b18434" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -17,6 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" DocStringExtensions = "~0.8, ~0.9" HDF5 = "0.17" InfrastructureSystems = "^1.20" +MKL = "0.6" KLU = "~0.4" Pardiso = "~0.5.4" PowerSystems = "3" diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl index abb3210d..1546cd00 100644 --- a/src/PowerNetworkMatrices.jl +++ b/src/PowerNetworkMatrices.jl @@ -19,6 +19,7 @@ export VirtualPTDF export Ybus using DocStringExtensions +import MKL import InfrastructureSystems import PowerSystems import PowerSystems: ACBusTypes @@ -31,9 +32,10 @@ import SparseArrays: rowvals, nzrange import HDF5 import KLU: klu import KLU -import LinearAlgebra: LAPACK.getri!, LAPACK.getrf!, BLAS.gemm -import LinearAlgebra: ldiv!, mul!, I, dot import LinearAlgebra +import LinearAlgebra: BLAS.gemm +import LinearAlgebra: ldiv!, mul!, I, dot +import LinearAlgebra: LAPACK.getrf!, LAPACK.getrs! import Pardiso @template DEFAULT = """ diff --git a/src/definitions.jl b/src/definitions.jl index 4941d0c0..b4d071a0 100644 --- a/src/definitions.jl +++ b/src/definitions.jl @@ -5,3 +5,5 @@ const MiB = KiB * KiB const GiB = MiB * KiB const MAX_CACHE_SIZE_MiB = 100 const ROW_PERSISTENT_CACHE_WARN = 1 * GiB + +DEFAULT_LODF_CHUNK_SIZE = 18_000 diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl index 79c5c061..98f12076 100644 --- a/src/lodf_calculations.jl +++ b/src/lodf_calculations.jl @@ -80,9 +80,9 @@ function _calculate_LODF_matrix_KLU( end end Dem_LU = klu(SparseArrays.sparse(m_I, m_I, m_V)) - lodf_t = Dem_LU \ ptdf_denominator - lodf_t[SparseArrays.diagind(lodf_t)] .= -1.0 - return lodf_t + KLU.solve!(Dem_LU, ptdf_denominator) + ptdf_denominator[SparseArrays.diagind(ptdf_denominator)] .= -1.0 + return ptdf_denominator end function _calculate_LODF_matrix_KLU( @@ -123,9 +123,101 @@ function _calculate_LODF_matrix_DENSE( push!(m_V, 1.0 - ptdf_denominator_t[iline, iline]) end end - lodf_t = LinearAlgebra.diagm(m_V) \ ptdf_denominator_t - lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0 - return lodf_t + (mV, bipiv, binfo) = getrf!(Matrix(LinearAlgebra.diagm(m_V))) + _binfo_check(binfo) + getrs!('N', mV, bipiv, ptdf_denominator_t) + ptdf_denominator_t[LinearAlgebra.diagind(ptdf_denominator_t)] .= -1.0 + return ptdf_denominator_t +end + +function _pardiso_sequential_LODF!( + lodf_t::Matrix{Float64}, + A::SparseArrays.SparseMatrixCSC{Float64, Int}, + ptdf_denominator_t::Matrix{Float64}, + chunk_size::Int = DEFAULT_LODF_CHUNK_SIZE, +) + @info "Line Count too large for single compute using Pardiso. Employing Sequential Calulations using a chunk_size=$(chunk_size)" + linecount = size(lodf_t, 1) + @assert LinearAlgebra.ishermitian(A) + ps = Pardiso.MKLPardisoSolver() + Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM) + Pardiso.pardisoinit(ps) + # Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON) + defaults = Pardiso.get_iparms(ps) + Pardiso.set_iparm!(ps, 1, 1) + for (ix, v) in enumerate(defaults[2:end]) + Pardiso.set_iparm!(ps, ix + 1, v) + end + Pardiso.set_iparm!(ps, 2, 2) + Pardiso.set_iparm!(ps, 59, 2) + Pardiso.set_iparm!(ps, 12, 1) + Pardiso.set_iparm!(ps, 11, 0) + Pardiso.set_iparm!(ps, 13, 0) + Pardiso.set_iparm!(ps, 32, 1) + #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON) + Pardiso.set_phase!(ps, Pardiso.ANALYSIS) + Pardiso.pardiso( + ps, + lodf_t, + A, + ptdf_denominator_t, + ) + + Pardiso.set_phase!(ps, Pardiso.NUM_FACT) + Pardiso.pardiso( + ps, + A, + Float64[], + ) + Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) + i_count = 1 + tmp = zeros(Float64, linecount, chunk_size) + while i_count <= linecount + edge = min(i_count + chunk_size - 1, linecount) + if linecount - edge <= 0 + tmp = tmp[:, 1:(edge - i_count + 1)] + end + Pardiso.pardiso( + ps, + tmp, + A, + ptdf_denominator_t[:, i_count:edge], + ) + lodf_t[:, i_count:edge] .= tmp + i_count = edge + 1 + end + Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) + Pardiso.pardiso(ps) + return +end + +function _pardiso_single_LODF!( + lodf_t::Matrix{Float64}, + A::SparseArrays.SparseMatrixCSC{Float64, Int}, + ptdf_denominator_t::Matrix{Float64}, +) + @assert LinearAlgebra.ishermitian(A) + ps = Pardiso.MKLPardisoSolver() + Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF) + Pardiso.pardisoinit(ps) + Pardiso.set_iparm!(ps, 1, 1) + defaults = Pardiso.get_iparms(ps) + for (ix, v) in enumerate(defaults[2:end]) + Pardiso.set_iparm!(ps, ix + 1, v) + end + Pardiso.set_iparm!(ps, 2, 2) + Pardiso.set_iparm!(ps, 59, 2) + Pardiso.set_iparm!(ps, 12, 1) + #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON) + Pardiso.pardiso( + ps, + lodf_t, + A, + ptdf_denominator_t, + ) + Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) + Pardiso.pardiso(ps) + return end function _calculate_LODF_matrix_MKLPardiso( @@ -145,28 +237,13 @@ function _calculate_LODF_matrix_MKLPardiso( push!(m_V, 1 - ptdf_denominator_t[iline, iline]) end end - - # intialize MKLPardiso - ps = Pardiso.MKLPardisoSolver() - defaults = Pardiso.get_iparms(ps) - Pardiso.set_iparm!(ps, 1, 1) - for (ix, v) in enumerate(defaults[2:end]) - Pardiso.set_iparm!(ps, ix + 1, v) - end - Pardiso.set_iparm!(ps, 2, 2) - Pardiso.set_iparm!(ps, 59, 2) - Pardiso.set_iparm!(ps, 6, 1) - # inizialize matrix for evaluation lodf_t = zeros(linecount, linecount) - # solve system - Pardiso.pardiso( - ps, - lodf_t, - SparseArrays.sparse(m_I, m_I, m_V), - ptdf_denominator_t, - ) - Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) - # set diagonal to -1 + A = SparseArrays.sparse(m_I, m_I, m_V) + if linecount > DEFAULT_LODF_CHUNK_SIZE + _pardiso_sequential_LODF!(lodf_t, A, ptdf_denominator_t) + else + _pardiso_single_LODF!(lodf_t, A, ptdf_denominator_t) + end lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0 return lodf_t end diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl index 558880b3..03e62a08 100644 --- a/src/ptdf_calculations.jl +++ b/src/ptdf_calculations.jl @@ -123,19 +123,17 @@ function _calculate_PTDF_matrix_KLU( # inizialize matrices for evaluation valid_ix = setdiff(1:buscount, ref_bus_positions) PTDFm_t = zeros(buscount, linecount) - + copyto!(PTDFm_t, BA) if !isempty(dist_slack) && length(ref_bus_positions) != 1 error( "Distibuted slack is not supported for systems with multiple reference buses.", ) elseif isempty(dist_slack) && length(ref_bus_positions) < buscount - copyto!(PTDFm_t, BA) PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :]) PTDFm_t[collect(ref_bus_positions), :] .= 0.0 return PTDFm_t elseif length(dist_slack) == buscount @info "Distributed bus" - copyto!(PTDFm_t, BA) PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :]) PTDFm_t[collect(ref_bus_positions), :] .= 0.0 slack_array = dist_slack / sum(dist_slack) @@ -172,6 +170,19 @@ function calculate_PTDF_matrix_KLU( return PTDFm, A end +function _binfo_check(binfo::Int) + if binfo != 0 + if binfo < 0 + error("Illegal Argument in Inputs") + elseif binfo > 0 + error("Singular value in factorization. Possibly there is an islanded bus") + else + @assert false + end + end + return +end + """ Funciton for internal use only. @@ -198,17 +209,21 @@ function _calculate_PTDF_matrix_DENSE( valid_ixs = setdiff(1:buscount, ref_bus_positions) ABA = Matrix(calculate_ABA_matrix(A, BA, ref_bus_positions)) PTDFm_t = zeros(buscount, linecount) - + (ABA, bipiv, binfo) = getrf!(ABA) + _binfo_check(binfo) + BA = Matrix(BA[valid_ixs, :]) if !isempty(dist_slack) && length(ref_bus_positions) != 1 error( "Distibuted slack is not supported for systems with multiple reference buses.", ) elseif isempty(dist_slack) && length(ref_bus_positions) < buscount - PTDFm_t[valid_ixs, :] = ABA \ BA[valid_ixs, :] + getrs!('N', ABA, bipiv, BA) + PTDFm_t[valid_ixs, :] = BA return PTDFm_t elseif length(dist_slack) == buscount @info "Distributed bus" - PTDFm_t[valid_ixs, :] = ABA \ BA[valid_ixs, :] + getrs!('N', ABA, bipiv, BA) + PTDFm_t[valid_ixs, :] = BA slack_array = dist_slack / sum(dist_slack) slack_array = reshape(slack_array, 1, buscount) return PTDFm_t - @@ -268,8 +283,10 @@ function _calculate_PTDF_matrix_MKLPardiso( buscount = size(BA, 1) ABA = calculate_ABA_matrix(A, BA, ref_bus_positions) - # Here add the subnetwork detection + @assert LinearAlgebra.issymmetric(ABA) ps = Pardiso.MKLPardisoSolver() + Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM) + Pardiso.pardisoinit(ps) # Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON) defaults = Pardiso.get_iparms(ps) Pardiso.set_iparm!(ps, 1, 1) @@ -279,6 +296,10 @@ function _calculate_PTDF_matrix_MKLPardiso( Pardiso.set_iparm!(ps, 2, 2) Pardiso.set_iparm!(ps, 59, 2) Pardiso.set_iparm!(ps, 6, 1) + Pardiso.set_iparm!(ps, 12, 1) + Pardiso.set_iparm!(ps, 11, 0) + Pardiso.set_iparm!(ps, 13, 0) + Pardiso.set_iparm!(ps, 32, 1) # inizialize matrices for evaluation valid_ix = setdiff(1:buscount, ref_bus_positions) @@ -293,12 +314,14 @@ function _calculate_PTDF_matrix_MKLPardiso( Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA) PTDFm_t[valid_ix, :] = full_BA Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) + Pardiso.pardiso(ps) return PTDFm_t elseif length(dist_slack) == buscount @info "Distributed bus" Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA) PTDFm_t[valid_ix, :] = full_BA Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) + Pardiso.pardiso(ps) slack_array = dist_slack / sum(dist_slack) slack_array = reshape(slack_array, 1, buscount) return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)