Merge pull request #59 from NREL-Sienna/afc/BLAS_improvements_and_MKL…

…Pardiso_LODF

Fixed Dense method and a new MKLPardiso for LODF evaluation

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
@@ -109,50 +109,65 @@ function _calculate_LODF_matrix_KLU(
     return lodf_t
 end
 
-# ! temp for evaluation
-function _calculate_LODF_matrix_KLU2(
+function _calculate_LODF_matrix_DENSE(
     a::SparseArrays.SparseMatrixCSC{Int8, Int},
     ptdf::Matrix{Float64},
 )
-    ptdf_denominator_t = a * ptdf
     linecount = size(ptdf, 2)
-    lodf_t = zeros(linecount, linecount)
-    for i in 1:linecount
-        for j in 1:linecount
-            if i == j
-                lodf_t[i, i] = -1.0
-            else
-                if (1.0 - ptdf_denominator_t[i, i]) < 1.0E-06
-                    lodf_t[i, j] = ptdf_denominator_t[i, j]
-                else
-                    lodf_t[i, j] = ptdf_denominator_t[i, j] / (1 - ptdf_denominator_t[i, i])
-                end
-            end
+    ptdf_denominator_t = a * ptdf
+    m_V = Float64[]
+    for iline in 1:linecount
+        if (1.0 - ptdf_denominator_t[iline, iline]) < 1.0E-06
+            push!(m_V, 1.0)
+        else
+            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
 end
 
-function _calculate_LODF_matrix_DENSE(
+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
-            ptdf_denominator_t[iline, iline] = 0.0
+            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
-    (Dem, dipiv, dinfo) = getrf!(
-        Matrix{Float64}(1.0I, linecount, linecount) -
-        Array(LinearAlgebra.Diagonal(ptdf_denominator_t)),
-    )
-    lodf_t = gemm('N', 'N', getri!(Dem, dipiv), ptdf_denominator_t)
-    lodf_t =
-        lodf_t - Array(LinearAlgebra.Diagonal(lodf_t)) -
-        Matrix{Float64}(1.0I, linecount, linecount)
 
+    # 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
 
@@ -164,6 +179,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())
 """

test/test_BA_ABA_matrix.jl

@@ -67,17 +67,7 @@ end
 
     # test if error is correctly thrown when ref bus is called
     rb = collect(ba.ref_bus_positions)[1]
-    test_val = false
-    try
-        aba[rb, :]
-    catch err
-        if err isa ErrorException
-            test_val = true
-        else
-            error("Expected an ErrorException but was not thrown.")
-        end
-    end
-    @test test_val
+    @test_throws ErrorException aba[rb, :]
 end
 
 @testset "Test show for A, BA and ABA matrix" begin

test/test_lodf.jl

@@ -13,12 +13,6 @@
     @test isapprox(maximum(L5.data), maximum(Lodf_5), atol = 1e-3)
     @test isapprox(L5[branches_5[1], branches_5[2]], 0.3447946513849093)
 
-    # get LODF with second method
-    a = IncidenceMatrix(sys5)
-    ptdf = PTDF(sys5)
-    lodf_t_3 = PNM._calculate_LODF_matrix_KLU2(a.data, ptdf.data)
-    @test isapprox(lodf_t_3, L5.data, atol = 1e-5)
-
     # get LODF with system
     L5NS = LODF(sys5)
     @test getindex(L5NS, "5", "6") - -0.3071 <= 1e-4
@@ -30,12 +24,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)
@@ -53,38 +51,14 @@
     end
 
     # test if error is thrown in case other linear solvers are called
-    test_value = false
-    try
-        L5NS_from_ba_aba = LODF(A, ABA, BA; linear_solver = "Dense")
-    catch err
-        if err isa ErrorException
-            test_value = true
-        end
-    end
-    @test test_value
+    @test_throws ErrorException LODF(A, ABA, BA; linear_solver = "Dense")
+
+    @test_throws ErrorException LODF(A, P5; linear_solver = "XXX")
 
     # 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)
-    catch err
-        if err isa ErrorException
-            test_value = true
-        end
-    end
-    @test test_value
-
-    # test if error is thrown in case `MKLPardiso` is chosen as a linera solver
-    test_value = false
-    try
-        lodf = LODF(sys5; linear_solver = "MKLPardiso")
-    catch err
-        if err isa ErrorException
-            test_value = true
-        end
-    end
-    @test test_value
+    @test_throws ErrorException L5NS_from_ptdf = LODF(A, P5)
 
     # get 14 bus system
     buses_14 = nodes14()

test/test_ptdf.jl

@@ -230,37 +230,13 @@ end
     dist_slack = 1 / buscount * ones(buscount)
     slack_array = dist_slack / sum(dist_slack)
 
-    test_val1 = false
-    try
-        ptdf_1 = PTDF(sys; dist_slack = slack_array)
-    catch err
-        if err isa ErrorException
-            test_val1 = true
-        else
-            error(
-                "Error was not thrown for PTDF with distributed slack bus and more than one reference bus.",
-            )
-        end
-    end
+    @test_throws ErrorException ptdf_1 = PTDF(sys; dist_slack = slack_array)
 
     # incorrect dist_slack arrya length
     sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5")
     buscount = length(PNM.get_buses(sys5)) + 1
     dist_slack = 1 / buscount * ones(buscount)
     slack_array = dist_slack / sum(dist_slack)
-    test_val2 = false
-    try
-        ptdf_2 = PTDF(sys5; dist_slack = slack_array)
-    catch err
-        if err isa ErrorException
-            test_val2 = true
-        else
-            error(
-                "Error was not thrown for PTDF with distributed slack array of incorrect length.",
-            )
-        end
-    end
 
-    @test test_val1
-    @test test_val2
+    @test_throws ErrorException ptdf_2 = PTDF(sys5; dist_slack = slack_array)
 end

test/test_virtual_lodf.jl

@@ -112,15 +112,7 @@ end
           length(axes(vlodf)[1]) * length(axes(vlodf)[2])
 
     # check if error is correctly thrown
-    test_value = false
-    try
-        vlodf[1, 1] = 1
-    catch err
-        if err isa ErrorException
-            test_value = true
-        end
-    end
-    @test test_value
+    @test_throws ErrorException vlodf[1, 1] = 1
 
     # test show
     test_value = false

test/test_virtual_ptdf.jl

@@ -101,41 +101,17 @@ end
     dist_slack = 1 / buscount * ones(buscount)
     slack_array = dist_slack / sum(dist_slack)
 
-    test_val1 = false
-    try
-        ptdf_1 = VirtualPTDF(sys; dist_slack = slack_array)
-        ptdf_1[1, 1]
-    catch err
-        if err isa ErrorException
-            test_val1 = true
-        else
-            error(
-                "Error was not thrown for PTDF with distributed slack bus and more than one reference bus.",
-            )
-        end
-    end
+    ptdf_1 = VirtualPTDF(sys; dist_slack = slack_array)
+    @test_throws ErrorException ptdf_1[1, 1]
 
     # incorrect dist_slack array length
     sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5")
     buscount = length(PNM.get_buses(sys5)) + 1
     dist_slack = 1 / buscount * ones(buscount)
     slack_array = dist_slack / sum(dist_slack)
     test_val2 = false
-    try
-        ptdf_2 = VirtualPTDF(sys5; dist_slack = slack_array)
-        ptdf_2[1, 1]
-    catch err
-        if err isa ErrorException
-            test_val2 = true
-        else
-            error(
-                "Error was not thrown for PTDF with distributed slack array of incorrect length.",
-            )
-        end
-    end
-
-    @test test_val1
-    @test test_val2
+    ptdf_2 = VirtualPTDF(sys5; dist_slack = slack_array)
+    @test_throws ErrorException ptdf_2[1, 1]
 end
 
 @testset "Test Virtual PTDF auxiliary functions" begin
@@ -150,15 +126,7 @@ end
     @test isempty(vptdf) == false
 
     # check if error is correctly thrown
-    test_value = false
-    try
-        vptdf[1, 1] = 1
-    catch err
-        if err isa ErrorException
-            test_value = true
-        end
-    end
-    @test test_value
+    @test_throws ErrorException vptdf[1, 1] = 1
 
     # get the rows and full PTDF matrix, test get_ptdf_data
     ptdf = PTDF(sys)