From b0409dc079dd75a6549cb5cac5a2bcd4573e8b62 Mon Sep 17 00:00:00 2001
From: alefcastelli <alef.castelli@gmail.com>
Date: Wed, 20 Sep 2023 22:12:50 -0600
Subject: [PATCH] LODF with MKLPardiso

---
 src/lodf_calculations.jl | 53 +++++++++++++++++++++++++++++++++++++---
 test/test_lodf.jl        | 13 ++++++----
 2 files changed, 58 insertions(+), 8 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index b69fb6bd..9a96a9ab 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -32,8 +32,8 @@ function _buildlodf(
         lodf_t = _calculate_LODF_matrix_KLU(a, ptdf)
     elseif linear_solver == "Dense"
         lodf_t = _calculate_LODF_matrix_DENSE(a, ptdf)
-    else
-        error("MKLPardiso still to be implemented")
+    elseif linear_solver == "MKLPardiso"
+        lodf_t = _calculate_LODF_matrix_MKLPardiso(a, ptdf)
     end
     return lodf_t
 end
@@ -147,7 +147,50 @@ function _calculate_LODF_matrix_DENSE(
             push!(m_V, 1.0 - ptdf_denominator_t[iline, iline])
         end
     end
-    lodf_t = LinearAlgebra.diagm(m_V) \ (a * ptdf)
+    lodf_t = LinearAlgebra.diagm(m_V) \ ptdf_denominator_t
+    lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
+    return lodf_t
+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
+
+    # 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
     lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
     return lodf_t
 end
@@ -160,6 +203,10 @@ Builds the LODF matrix given the data of branches and buses of the system.
         vector of the System AC branches
 - `buses::Vector{PSY.ACBus}`:
         vector of the System buses
+- `linear_solver`::String
+        linear solver to use for matrix evaluation.
+        Available options: "KLU", "Dense" (OpenBLAS) or "MKLPardiso".
+        Default solver: "KLU".
 - `tol::Float64`:
         Tolerance to eliminate entries in the LODF matrix (default eps())
 """
diff --git a/test/test_lodf.jl b/test/test_lodf.jl
index 5a26cc2a..8b6bea15 100644
--- a/test/test_lodf.jl
+++ b/test/test_lodf.jl
@@ -30,12 +30,16 @@
     P5 = PTDF(sys5)
     L5NS_from_ptdf = LODF(A, P5)
     L5NS_from_ptdf2 = LODF(A, P5; linear_solver = "Dense")
+    L5NS_from_ptdf3 = LODF(A, P5; linear_solver = "MKLPardiso")
     @test getindex(L5NS_from_ptdf, "5", "6") - -0.3071 <= 1e-4
     @test getindex(L5NS_from_ptdf2, "5", "6") - -0.3071 <= 1e-4
+    @test getindex(L5NS_from_ptdf3, "5", "6") - -0.3071 <= 1e-4
     total_error = abs.(L5NS_from_ptdf.data' .- Lodf_5)
     total_error2 = abs.(L5NS_from_ptdf2.data' .- Lodf_5)
+    total_error3 = abs.(L5NS_from_ptdf3.data' .- Lodf_5)
     @test isapprox(sum(total_error), 0.0, atol = 1e-3)
     @test isapprox(sum(total_error2), 0.0, atol = 1e-3)
+    @test isapprox(sum(total_error3), 0.0, atol = 1e-3)
 
     # A, ABA, and BA case
     ABA = ABA_Matrix(sys5; factorize = true)
@@ -63,11 +67,9 @@
     end
     @test test_value
 
-    # test if error is thrown in case `tol` is defined in PTDF
-    P5 = PTDF(sys5; tol = 1e-3)
     test_value = false
     try
-        L5NS_from_ptdf = LODF(A, P5)
+        L5NS_from_ba_aba = LODF(A, P5; linear_solver = "XXX")
     catch err
         if err isa ErrorException
             test_value = true
@@ -75,10 +77,11 @@
     end
     @test test_value
 
-    # test if error is thrown in case `MKLPardiso` is chosen as a linera solver
+    # test if error is thrown in case `tol` is defined in PTDF
+    P5 = PTDF(sys5; tol = 1e-3)
     test_value = false
     try
-        lodf = LODF(sys5; linear_solver = "MKLPardiso")
+        L5NS_from_ptdf = LODF(A, P5)
     catch err
         if err isa ErrorException
             test_value = true