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

src/lodf_calculations.jl

@@ -128,6 +128,73 @@
     return lodf_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,
+)
+    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
+
+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()
+    Pardiso.pardisoinit(ps)
+    #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    Pardiso.pardiso(
+        ps,
+        lodf_t,
+        A,
+        ptdf_denominator_t,
+    )
+    Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+    return
+end
+
+
 function _calculate_LODF_matrix_MKLPardiso(
     a::SparseArrays.SparseMatrixCSC{Int8, Int},
     ptdf::Matrix{Float64},
@@ -145,28 +212,13 @@
             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