Merge pull request #56 from NREL-Sienna/jd/debug_pardiso

Jd/debug pardiso

Project.toml

@@ -18,6 +18,6 @@ DocStringExtensions = "~0.8, ~0.9"
 HDF5 = "0.17"
 InfrastructureSystems = "^1.20"
 KLU = "~0.4"
-Pardiso = "~0.5"
+Pardiso = "~0.5.4"
 PowerSystems = "3"
 julia = "^1.6"

src/ptdf_calculations.jl

@@ -295,38 +295,42 @@ function _calculate_PTDF_matrix_MKLPardiso(
 
     ABA = calculate_ABA_matrix(A, BA, ref_bus_positions)
     # Here add the subnetwork detection
-    Ix = Matrix(
-        1.0I,
-        buscount - length(ref_bus_positions),
-        buscount - length(ref_bus_positions),
-    )
-    ABA_inv = zeros(Float64, size(Ix))
-
     ps = Pardiso.MKLPardisoSolver()
+    # 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_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
-    Pardiso.solve!(ps, ABA_inv, ABA, Ix)
+    Pardiso.set_iparm!(ps, 6, 1)
+
     # inizialize matrices for evaluation
     valid_ix = setdiff(1:buscount, ref_bus_positions)
     PTDFm_t = zeros(buscount, linecount)
 
+    full_BA = Matrix(BA[valid_ix, :])
     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_ix, :] = ABA_inv * @view BA[valid_ix, :]
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        PTDFm_t[valid_ix, :] = full_BA
+        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        PTDFm_t[valid_ix, :] = ABA_inv * @view BA[valid_ix, :]
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        PTDFm_t[valid_ix, :] = full_BA
+        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)
     else
         error("Distributed bus specification doesn't match the number of buses.")
     end
-
     return
 end
 

test/Project.toml

@@ -6,6 +6,7 @@ InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"

test/test_ptdf.jl

@@ -200,7 +200,7 @@ end
 
 @testset "Test PTDF matrices with distributed slack" begin
     """
-    CAUTION: this test just test that all the matrices are the same, but there 
+    CAUTION: this test just test that all the matrices are the same, but there
     are no reference values.
     """
 
@@ -215,9 +215,9 @@ end
     P5_2 = PTDF(sys5; dist_slack = slack_array, linear_solver = "Dense")
     P5_3 = PTDF(sys5; dist_slack = slack_array, linear_solver = "MKLPardiso")
 
-    @assert isapprox(P5_1.data, P5_2.data, atol = 1e-5)
-    @assert isapprox(P5_1.data, P5_3.data, atol = 1e-5)
-    @assert isapprox(P5_2.data, P5_3.data, atol = 1e-5)
+    @test isapprox(P5_1.data, P5_2.data, atol = 1e-5)
+    @test isapprox(P5_1.data, P5_3.data, atol = 1e-5)
+    @test isapprox(P5_2.data, P5_3.data, atol = 1e-5)
 end
 
 @testset "Test PTDF matrix with distributed bus and with 2 reference buses" begin