Jd/debug lodf

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 = 20_000

src/lodf_calculations.jl

@@ -128,45 +128,105 @@ function _calculate_LODF_matrix_DENSE(
     return lodf_t
 end
 
-function _calculate_LODF_matrix_MKLPardiso(
-    a::SparseArrays.SparseMatrixCSC{Int8, Int},
-    ptdf::Matrix{Float64},
+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,
 )
-    linecount = size(ptdf, 2)
-    ptdf_denominator_t = a * ptdf
-    m_I = Int[]
-    m_V = Float64[]
-    for iline in 1:linecount
-        if (1.0 - ptdf_denominator_t[iline, iline]) < 1.0E-06
-            push!(m_I, iline)
-            push!(m_V, 1.0)
-        else
-            push!(m_I, iline)
-            push!(m_V, 1 - ptdf_denominator_t[iline, iline])
+    linecount = size(lodf_t, 1)
+    ps = Pardiso.MKLPardisoSolver()
+    Pardiso.pardisoinit(ps)
+    #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)
+    return
+end
 
-    # intialize MKLPardiso
+function _pardiso_single_LODF!(
+    lodf_t::Matrix{Float64},
+    A::SparseArrays.SparseMatrixCSC{Float64, Int},
+    ptdf_denominator_t::Matrix{Float64},
+    chunk_size::Int = DEFAULT_LODF_CHUNK_SIZE,
+)
+    linecount = size(lodf_t, 1)
     ps = Pardiso.MKLPardisoSolver()
-    defaults = Pardiso.get_iparms(ps)
+    Pardiso.pardisoinit(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.set_iparm!(ps, 12, 1)
+    #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     Pardiso.pardiso(
         ps,
         lodf_t,
-        SparseArrays.sparse(m_I, m_I, m_V),
+        A,
         ptdf_denominator_t,
     )
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-    # set diagonal to -1
+    return
+end
+
+
+function _calculate_LODF_matrix_MKLPardiso(
+    a::SparseArrays.SparseMatrixCSC{Int8, Int},
+    ptdf::Matrix{Float64},
+)
+    linecount = size(ptdf, 2)
+    ptdf_denominator_t = a * ptdf
+    m_I = Int[]
+    m_V = Float64[]
+    for iline in 1:linecount
+        if (1.0 - ptdf_denominator_t[iline, iline]) < 1.0E-06
+            push!(m_I, iline)
+            push!(m_V, 1.0)
+        else
+            push!(m_I, iline)
+            push!(m_V, 1 - ptdf_denominator_t[iline, iline])
+        end
+    end
+    lodf_t = zeros(linecount, linecount)
+    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

@@ -270,6 +270,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
     ABA = calculate_ABA_matrix(A, BA, ref_bus_positions)
     # Here add the subnetwork detection
     ps = Pardiso.MKLPardisoSolver()
+    Pardiso.pardisoinit(ps)
     # Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     defaults = Pardiso.get_iparms(ps)
     Pardiso.set_iparm!(ps, 1, 1)
@@ -279,6 +280,7 @@ 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)
 
     # inizialize matrices for evaluation
     valid_ix = setdiff(1:buscount, ref_bus_positions)
@@ -290,15 +292,17 @@ function _calculate_PTDF_matrix_MKLPardiso(
             "Distibuted slack is not supported for systems with multiple reference buses.",
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
-        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
-        PTDFm_t[valid_ix, :] = full_BA
+        Pardiso.pardiso(ps, similar(full_BA), ABA, full_BA)
+        PTDFm_t[valid_ix, :] .= full_BA
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.pardiso(ps, Float64[], ABA, full_BA)
         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.pardiso(ps, similar(full_BA), ABA, full_BA)
+        PTDFm_t[valid_ix, :] .= full_BA
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)