Merge pull request #58 from NREL-Sienna/BLAS_improvements

modified PTDF functions for Dense method

Project.toml

@@ -1,6 +1,6 @@
 name = "PowerNetworkMatrices"
 uuid = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
-authors = ["Jose Daniel Lara", "Alessandro Castelli", "Sourabh Dalvi"]
+authors = ["Jose Daniel Lara", "Alessandro Francesco Castelli", "Sourabh Dalvi"]
 version = "0.9.3"
 
 [deps]

src/ptdf_calculations.jl

@@ -172,19 +172,6 @@ function calculate_PTDF_matrix_KLU(
     return PTDFm, A
 end
 
-function _binfo_check(binfo::Int)
-    if binfo != 0
-        if binfo < 0
-            error("Illegal Argument in Inputs")
-        elseif binfo > 0
-            error("Singular value in factorization. Possibly there is an islanded bus")
-        else
-            @assert false
-        end
-    end
-    return
-end
-
 """
 Funciton for internal use only.
 
@@ -201,40 +188,27 @@ Computes the PTDF matrix by means of the LAPACK and BLAS functions for dense mat
         vector containing the weights for the distributed slacks.
 """
 function _calculate_PTDF_matrix_DENSE(
-    A::Matrix{Int8},
-    BA::Matrix{T},
+    A::SparseArrays.SparseMatrixCSC{Int8, Int},
+    BA::SparseArrays.SparseMatrixCSC{Float64, Int},
     ref_bus_positions::Set{Int},
-    dist_slack::Vector{Float64}) where {T <: Union{Float32, Float64}}
+    dist_slack::Vector{Float64})
     linecount = size(BA, 2)
     buscount = size(BA, 1)
     # Use dense calculation of ABA
     valid_ixs = setdiff(1:buscount, ref_bus_positions)
-    ABA = transpose(A[:, valid_ixs]) * transpose(BA[valid_ixs, :])
-    # get LU factorization matrices
-    (ABA, bipiv, binfo) = getrf!(ABA)
-    _binfo_check(binfo)
+    ABA = Matrix(calculate_ABA_matrix(A, BA, ref_bus_positions))
     PTDFm_t = zeros(buscount, linecount)
 
     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_ixs, :] = gemm(
-            'N',
-            'N',
-            getri!(ABA, bipiv),
-            BA[valid_ixs, :],
-        )
+        PTDFm_t[valid_ixs, :] = ABA \ BA[valid_ixs, :]
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        PTDFm_t[valid_ixs, :] = gemm(
-            'N',
-            'N',
-            getri!(ABA, bipiv),
-            BA[valid_ixs, :],
-        )
+        PTDFm_t[valid_ixs, :] = ABA \ BA[valid_ixs, :]
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t -
@@ -265,8 +239,8 @@ function calculate_PTDF_matrix_DENSE(
     bus_lookup::Dict{Int, Int},
     dist_slack::Vector{Float64})
     A, ref_bus_positions = calculate_A_matrix(branches, buses)
-    BA = Matrix(calculate_BA_matrix(branches, bus_lookup))
-    PTDFm = _calculate_PTDF_matrix_DENSE(Matrix(A), BA, ref_bus_positions, dist_slack)
+    BA = calculate_BA_matrix(branches, bus_lookup)
+    PTDFm = _calculate_PTDF_matrix_DENSE(A, BA, ref_bus_positions, dist_slack)
     return PTDFm, A
 end