Jd/debug lodf

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"

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 = """

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

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

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)