From 51abb8a20ca1022ae5279a617ad3e42293f1b028 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 12:13:24 -0600
Subject: [PATCH 01/56] add debugging code

---
 src/lodf_calculations.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 79c5c061..d1f1a287 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -132,6 +132,7 @@ function _calculate_LODF_matrix_MKLPardiso(
     a::SparseArrays.SparseMatrixCSC{Int8, Int},
     ptdf::Matrix{Float64},
 )
+    @error "LODF Start"
     linecount = size(ptdf, 2)
     ptdf_denominator_t = a * ptdf
     m_I = Int[]
@@ -159,6 +160,8 @@ function _calculate_LODF_matrix_MKLPardiso(
     # inizialize matrix for evaluation
     lodf_t = zeros(linecount, linecount)
     # solve system
+    Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    @error "Call to Pardiso Start"
     Pardiso.pardiso(
         ps,
         lodf_t,

From ea07a25d5ce6c903a14ed9027f805c65b2925fd2 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 13:51:23 -0600
Subject: [PATCH 02/56] add more checking

---
 src/lodf_calculations.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index d1f1a287..ab201a4a 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -157,6 +157,7 @@ function _calculate_LODF_matrix_MKLPardiso(
     Pardiso.set_iparm!(ps, 2, 2)
     Pardiso.set_iparm!(ps, 59, 2)
     Pardiso.set_iparm!(ps, 6, 1)
+    Pardiso.set_iparm!(ps, 27, 1)
     # inizialize matrix for evaluation
     lodf_t = zeros(linecount, linecount)
     # solve system

From 6465028ef61868d64de1849291abb68d6f7054a0 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 15:51:50 -0600
Subject: [PATCH 03/56] one more test

---
 src/lodf_calculations.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index ab201a4a..690c8748 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -151,13 +151,16 @@ function _calculate_LODF_matrix_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
+    #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)
     Pardiso.set_iparm!(ps, 27, 1)
+    Pardiso.set_iparm!(ps, 10, 0)
+    Pardiso.set_iparm!(ps, 12, 0)
+    Pardiso.set_iparm!(ps, 23, 0)
     # inizialize matrix for evaluation
     lodf_t = zeros(linecount, linecount)
     # solve system

From 6556db888b0b956df16252f7acc1648fd25ed29e Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 16:07:17 -0600
Subject: [PATCH 04/56] remove write to vector option

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 690c8748..1edd3e1d 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -156,7 +156,7 @@ function _calculate_LODF_matrix_MKLPardiso(
     #end
     Pardiso.set_iparm!(ps, 2, 2)
     Pardiso.set_iparm!(ps, 59, 2)
-    Pardiso.set_iparm!(ps, 6, 1)
+    # Pardiso.set_iparm!(ps, 6, 1)
     Pardiso.set_iparm!(ps, 27, 1)
     Pardiso.set_iparm!(ps, 10, 0)
     Pardiso.set_iparm!(ps, 12, 0)

From 2fcb7a4047f54d71fe567b62eccaa6406d51716a Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 16:26:45 -0600
Subject: [PATCH 05/56] split phases

---
 src/lodf_calculations.jl | 22 ++++++++++++++++++++--
 1 file changed, 20 insertions(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 1edd3e1d..92976c09 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -149,11 +149,14 @@ function _calculate_LODF_matrix_MKLPardiso(
 
     # intialize MKLPardiso
     ps = Pardiso.MKLPardisoSolver()
+    Pardiso.pardisoinit(ps)
+    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, 12, 1)
     Pardiso.set_iparm!(ps, 2, 2)
     Pardiso.set_iparm!(ps, 59, 2)
     # Pardiso.set_iparm!(ps, 6, 1)
@@ -164,8 +167,23 @@ function _calculate_LODF_matrix_MKLPardiso(
     # inizialize matrix for evaluation
     lodf_t = zeros(linecount, linecount)
     # solve system
-    Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
-    @error "Call to Pardiso Start"
+    @error "Call to Pardiso Analysis"
+    Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
+    Pardiso.pardiso(
+        ps,
+        lodf_t,
+        SparseArrays.sparse(m_I, m_I, m_V),
+        ptdf_denominator_t,
+    )
+    Pardiso.set_phase!(ps, Pardiso.NUM_FACT)
+    @error "Call to Pardiso Fact"
+    Pardiso.pardiso(
+        ps,
+        ptdf_denominator_t,
+        Float64[]
+    )
+    Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE)
+    @error "Call to Pardiso Solve"
     Pardiso.pardiso(
         ps,
         lodf_t,

From 3fe37e706451a0ecd23a169a021df018d23248ea Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 16:43:19 -0600
Subject: [PATCH 06/56] fix typo

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 92976c09..6287332d 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -179,7 +179,7 @@ function _calculate_LODF_matrix_MKLPardiso(
     @error "Call to Pardiso Fact"
     Pardiso.pardiso(
         ps,
-        ptdf_denominator_t,
+        SparseArrays.sparse(m_I, m_I, m_V),
         Float64[]
     )
     Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE)

From 5904898b71823a76a9d4b2036a958f1ad8f522d8 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 16:54:12 -0600
Subject: [PATCH 07/56] fix typo

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 6287332d..a1aebb8d 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -182,7 +182,7 @@ function _calculate_LODF_matrix_MKLPardiso(
         SparseArrays.sparse(m_I, m_I, m_V),
         Float64[]
     )
-    Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE)
+    Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
     Pardiso.pardiso(
         ps,

From 4825d8b513460b6c394fe96e1f77143ff81ed4e3 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 19:15:02 -0600
Subject: [PATCH 08/56] try single solve call

---
 src/lodf_calculations.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index a1aebb8d..ca026f69 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -188,9 +188,10 @@ function _calculate_LODF_matrix_MKLPardiso(
         ps,
         lodf_t,
         SparseArrays.sparse(m_I, m_I, m_V),
-        ptdf_denominator_t,
+        ptdf_denominator_t[:,1],
     )
-    Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+    error()
+    # Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
     # set diagonal to -1
     lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
     return lodf_t

From cdbed5d1425eacbbed3c7584adc09bf021d86a7b Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 22:09:58 -0600
Subject: [PATCH 09/56] fix index

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index ca026f69..9f607f35 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -186,7 +186,7 @@ function _calculate_LODF_matrix_MKLPardiso(
     @error "Call to Pardiso Solve"
     Pardiso.pardiso(
         ps,
-        lodf_t,
+        lodf_t[:,1],
         SparseArrays.sparse(m_I, m_I, m_V),
         ptdf_denominator_t[:,1],
     )

From b40b98bcde5c0f9709ae0821df093d0ebb66d754 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 22:34:00 -0600
Subject: [PATCH 10/56] try total

---
 src/lodf_calculations.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 9f607f35..ba6ef1fb 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -186,11 +186,10 @@ function _calculate_LODF_matrix_MKLPardiso(
     @error "Call to Pardiso Solve"
     Pardiso.pardiso(
         ps,
-        lodf_t[:,1],
+        lodf_t,
         SparseArrays.sparse(m_I, m_I, m_V),
-        ptdf_denominator_t[:,1],
+        ptdf_denominator_t,
     )
-    error()
     # Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
     # set diagonal to -1
     lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0

From 9d8fe72a9581b4b2e5f24f5b8a3210ca1954eaab Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 22:57:04 -0600
Subject: [PATCH 11/56] test sequential solves

---
 src/lodf_calculations.jl | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index ba6ef1fb..ab055a0b 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -184,15 +184,20 @@ function _calculate_LODF_matrix_MKLPardiso(
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
-    Pardiso.pardiso(
-        ps,
-        lodf_t,
-        SparseArrays.sparse(m_I, m_I, m_V),
-        ptdf_denominator_t,
-    )
+    for i in 1:linecount
+        temp = zeros(Float64, linecount)
+        Pardiso.pardiso(
+            ps,
+            temp,
+            SparseArrays.sparse(m_I, m_I, m_V),
+            ptdf_denominator_t[:, i],
+        )
+        lodf_t[:,i] .= temp
+        lodf_t[i,i] = -1.0
+    end
     # Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
     # set diagonal to -1
-    lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
+    #lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
     return lodf_t
 end
 

From 01e4496e4af4df3b8ced8d743af3e93937365355 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 23:28:03 -0600
Subject: [PATCH 12/56] next try

---
 src/lodf_calculations.jl | 33 +++++++++++++++------------------
 1 file changed, 15 insertions(+), 18 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index ab055a0b..c3ed9f59 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -157,47 +157,44 @@ function _calculate_LODF_matrix_MKLPardiso(
     #    Pardiso.set_iparm!(ps, ix + 1, v)
     #end
     Pardiso.set_iparm!(ps, 12, 1)
-    Pardiso.set_iparm!(ps, 2, 2)
+    Pardiso.set_iparm!(ps, 2, 3)
     Pardiso.set_iparm!(ps, 59, 2)
     # Pardiso.set_iparm!(ps, 6, 1)
     Pardiso.set_iparm!(ps, 27, 1)
-    Pardiso.set_iparm!(ps, 10, 0)
-    Pardiso.set_iparm!(ps, 12, 0)
+    Pardiso.set_iparm!(ps, 11, 1)
+    Pardiso.set_iparm!(ps, 13, 1)
     Pardiso.set_iparm!(ps, 23, 0)
     # inizialize matrix for evaluation
     lodf_t = zeros(linecount, linecount)
     # solve system
     @error "Call to Pardiso Analysis"
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
+    A = SparseArrays.sparse(m_I, m_I, m_V)
     Pardiso.pardiso(
         ps,
         lodf_t,
-        SparseArrays.sparse(m_I, m_I, m_V),
+        A,
         ptdf_denominator_t,
     )
     Pardiso.set_phase!(ps, Pardiso.NUM_FACT)
     @error "Call to Pardiso Fact"
     Pardiso.pardiso(
         ps,
-        SparseArrays.sparse(m_I, m_I, m_V),
+        A,
         Float64[]
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
-    for i in 1:linecount
-        temp = zeros(Float64, linecount)
-        Pardiso.pardiso(
-            ps,
-            temp,
-            SparseArrays.sparse(m_I, m_I, m_V),
-            ptdf_denominator_t[:, i],
-        )
-        lodf_t[:,i] .= temp
-        lodf_t[i,i] = -1.0
-    end
-    # Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+    Pardiso.pardiso(
+        ps,
+        lodf_t,
+        A,
+        lodf_t,
+    )
+    @error "Call to Pardiso release"
+    Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
     # set diagonal to -1
-    #lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
+    lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
     return lodf_t
 end
 

From 550fdd3afb784596e376fb91d7f97171287d02e6 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Tue, 26 Sep 2023 23:48:47 -0600
Subject: [PATCH 13/56] use correct stage

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index c3ed9f59..53e49758 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -189,7 +189,7 @@ function _calculate_LODF_matrix_MKLPardiso(
         ps,
         lodf_t,
         A,
-        lodf_t,
+        ptdf_denominator_t,
     )
     @error "Call to Pardiso release"
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)

From b2e08c1148f7f281cadc5a0555f27afa06b5558a Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 00:18:04 -0600
Subject: [PATCH 14/56] test using chunks

---
 src/lodf_calculations.jl | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 53e49758..e1ad415c 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -185,12 +185,18 @@ function _calculate_LODF_matrix_MKLPardiso(
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
-    Pardiso.pardiso(
-        ps,
-        lodf_t,
-        A,
-        ptdf_denominator_t,
-    )
+    i_count = 1
+    while i_count <= linecount
+        edge = min(i_count + 2999, linecount)
+        tmp = zeros(Float64, 3000)
+        Pardiso.pardiso(
+            ps,
+            tmp,
+            A,
+            ptdf_denominator_t[:, i_count:edge],
+        )
+        i_count = edge
+    end
     @error "Call to Pardiso release"
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
     # set diagonal to -1

From 92dd80f82e6bd4f94398f793d7c306e9a5b5b3d6 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 00:27:43 -0600
Subject: [PATCH 15/56] update calculations

---
 src/lodf_calculations.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index e1ad415c..d70ea587 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -186,15 +186,16 @@ function _calculate_LODF_matrix_MKLPardiso(
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
     i_count = 1
+    tmp = zeros(Float64, 3000)
     while i_count <= linecount
         edge = min(i_count + 2999, linecount)
-        tmp = zeros(Float64, 3000)
         Pardiso.pardiso(
             ps,
             tmp,
             A,
             ptdf_denominator_t[:, i_count:edge],
         )
+        lodf_t[:, i_count:edge] .= tmp
         i_count = edge
     end
     @error "Call to Pardiso release"

From 6c35043858892aefd82c14f33a282c52837b1710 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 00:28:23 -0600
Subject: [PATCH 16/56] increase chunk size

---
 src/lodf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index d70ea587..f6f2e9fd 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -186,9 +186,9 @@ function _calculate_LODF_matrix_MKLPardiso(
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
     i_count = 1
-    tmp = zeros(Float64, 3000)
+    tmp = zeros(Float64, 10000)
     while i_count <= linecount
-        edge = min(i_count + 2999, linecount)
+        edge = min(i_count + 9999, linecount)
         Pardiso.pardiso(
             ps,
             tmp,

From 93cc28b2e470c09a4d3179b121f97027b706c87c Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 00:34:53 -0600
Subject: [PATCH 17/56] fix indexing

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index f6f2e9fd..62fff844 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -196,7 +196,7 @@ function _calculate_LODF_matrix_MKLPardiso(
             ptdf_denominator_t[:, i_count:edge],
         )
         lodf_t[:, i_count:edge] .= tmp
-        i_count = edge
+        i_count = edge + 1
     end
     @error "Call to Pardiso release"
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)

From 220c25b140a38ef8d9b2f5653397a6d27aeafda8 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 00:45:44 -0600
Subject: [PATCH 18/56] fix tmp size

---
 src/lodf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 62fff844..da8fd257 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -186,7 +186,7 @@ function _calculate_LODF_matrix_MKLPardiso(
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     @error "Call to Pardiso Solve"
     i_count = 1
-    tmp = zeros(Float64, 10000)
+    tmp = zeros(Float64, linecount, 10000)
     while i_count <= linecount
         edge = min(i_count + 9999, linecount)
         Pardiso.pardiso(

From de38f02fb50bcb5915c72acd11bbfc41d23124e4 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 01:04:50 -0600
Subject: [PATCH 19/56] address last chunk

---
 src/lodf_calculations.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index da8fd257..a19d6cb9 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -189,6 +189,9 @@ function _calculate_LODF_matrix_MKLPardiso(
     tmp = zeros(Float64, linecount, 10000)
     while i_count <= linecount
         edge = min(i_count + 9999, linecount)
+        if edge - i_count < 10000
+            tmp = tmp[:, i_count:edge]
+        end
         Pardiso.pardiso(
             ps,
             tmp,

From 35ca5753672cfb42e9d2ec1cfb61764742819f67 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 14:30:40 -0600
Subject: [PATCH 20/56] add sequential Pardiso solve

---
 src/definitions.jl       |  2 +
 src/lodf_calculations.jl | 84 ++++++++++++++++++----------------------
 2 files changed, 40 insertions(+), 46 deletions(-)

diff --git a/src/definitions.jl b/src/definitions.jl
index 4941d0c0..1501c6ff 100644
--- a/src/definitions.jl
+++ b/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
diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index a19d6cb9..786d92cc 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -128,69 +128,39 @@ function _calculate_LODF_matrix_DENSE(
     return lodf_t
 end
 
-function _calculate_LODF_matrix_MKLPardiso(
-    a::SparseArrays.SparseMatrixCSC{Int8, Int},
-    ptdf::Matrix{Float64},
-)
-    @error "LODF Start"
-    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
+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)
-    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, 12, 1)
-    Pardiso.set_iparm!(ps, 2, 3)
-    Pardiso.set_iparm!(ps, 59, 2)
-    # Pardiso.set_iparm!(ps, 6, 1)
-    Pardiso.set_iparm!(ps, 27, 1)
-    Pardiso.set_iparm!(ps, 11, 1)
-    Pardiso.set_iparm!(ps, 13, 1)
-    Pardiso.set_iparm!(ps, 23, 0)
-    # inizialize matrix for evaluation
-    lodf_t = zeros(linecount, linecount)
-    # solve system
-    @error "Call to Pardiso Analysis"
+
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
-    A = SparseArrays.sparse(m_I, m_I, m_V)
     Pardiso.pardiso(
         ps,
         lodf_t,
         A,
         ptdf_denominator_t,
     )
+
     Pardiso.set_phase!(ps, Pardiso.NUM_FACT)
-    @error "Call to Pardiso Fact"
     Pardiso.pardiso(
         ps,
         A,
         Float64[]
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
-    @error "Call to Pardiso Solve"
     i_count = 1
-    tmp = zeros(Float64, linecount, 10000)
+    tmp = zeros(Float64, linecount, chunk_size)
+
     while i_count <= linecount
-        edge = min(i_count + 9999, linecount)
-        if edge - i_count < 10000
-            tmp = tmp[:, i_count:edge]
+        edge = min(i_count + chunk_size, linecount)
+        if linecount - edge <= 0
+            tmp = tmp[:, 1:(edge - i_count + 1)]
         end
         Pardiso.pardiso(
             ps,
@@ -201,9 +171,31 @@ function _calculate_LODF_matrix_MKLPardiso(
         lodf_t[:, i_count:edge] .= tmp
         i_count = edge + 1
     end
-    @error "Call to Pardiso release"
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-    # set diagonal to -1
+    return
+end
+
+function _calculate_LODF_matrix_MKLPardiso(
+    a::SparseArrays.SparseMatrixCSC{Int8, Int},
+    ptdf::Matrix{Float64},
+)
+    @error "LODF Start"
+    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
+    lodf_t = zeros(linecount, linecount)
+    A = SparseArrays.sparse(m_I, m_I, m_V)
+    _pardiso_sequential_LODF!(lodf_t, A, ptdf_denominator_t)
     lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
     return lodf_t
 end

From 68f0e7621f86dc212bfd10533c1191a485c47951 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 14:32:56 -0600
Subject: [PATCH 21/56] formatter

---
 src/lodf_calculations.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 786d92cc..c9556c2e 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -132,8 +132,8 @@ 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
-    )
+    chunk_size::Int = DEFAULT_LODF_CHUNK_SIZE,
+)
     linecount = size(lodf_t, 1)
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.pardisoinit(ps)
@@ -151,7 +151,7 @@ function _pardiso_sequential_LODF!(
     Pardiso.pardiso(
         ps,
         A,
-        Float64[]
+        Float64[],
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
     i_count = 1

From 7b6a8c77c065f8fa48cf31d043620ceebf7bf22e Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 14:51:01 -0600
Subject: [PATCH 22/56] fix indexing

---
 src/lodf_calculations.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index c9556c2e..dd78a4ab 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -156,9 +156,8 @@ function _pardiso_sequential_LODF!(
     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, linecount)
+        edge = min(i_count + chunk_size - 1, linecount)
         if linecount - edge <= 0
             tmp = tmp[:, 1:(edge - i_count + 1)]
         end

From cde03672468cdf1e54fbca8178bab5ffc0f84c7b Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 14:51:21 -0600
Subject: [PATCH 23/56] count on 0

---
 src/lodf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index dd78a4ab..a1932bdf 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -154,10 +154,10 @@ function _pardiso_sequential_LODF!(
         Float64[],
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
-    i_count = 1
+    i_count = 0
     tmp = zeros(Float64, linecount, chunk_size)
     while i_count <= linecount
-        edge = min(i_count + chunk_size - 1, linecount)
+        edge = min(i_count + chunk_size, linecount)
         if linecount - edge <= 0
             tmp = tmp[:, 1:(edge - i_count + 1)]
         end

From 2886e341f2d1223425faa36fad7a8fafa43f3039 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 14:59:15 -0600
Subject: [PATCH 24/56] one more fix to indexing

---
 src/lodf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index a1932bdf..55e55f35 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -154,10 +154,10 @@ function _pardiso_sequential_LODF!(
         Float64[],
     )
     Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
-    i_count = 0
+    i_count = 1
     tmp = zeros(Float64, linecount, chunk_size)
     while i_count <= linecount
-        edge = min(i_count + chunk_size, linecount)
+        edge = min(i_count + chunk_size -  1, linecount)
         if linecount - edge <= 0
             tmp = tmp[:, 1:(edge - i_count + 1)]
         end

From 5f7736c66d05067b32770d009b55a26608603c65 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 18:58:35 -0600
Subject: [PATCH 25/56] remove logging

---
 src/lodf_calculations.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 55e55f35..fba94de1 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -137,7 +137,7 @@ function _pardiso_sequential_LODF!(
     linecount = size(lodf_t, 1)
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.pardisoinit(ps)
-    Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
 
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
     Pardiso.pardiso(
@@ -178,7 +178,6 @@ function _calculate_LODF_matrix_MKLPardiso(
     a::SparseArrays.SparseMatrixCSC{Int8, Int},
     ptdf::Matrix{Float64},
 )
-    @error "LODF Start"
     linecount = size(ptdf, 2)
     ptdf_denominator_t = a * ptdf
     m_I = Int[]

From 0d78a5d810ecde2f85941c884b60628e04f0a89f Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 19:02:00 -0600
Subject: [PATCH 26/56] enable single call

---
 src/lodf_calculations.jl | 27 ++++++++++++++++++++++++++-
 1 file changed, 26 insertions(+), 1 deletion(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index fba94de1..cf849d0b 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -174,6 +174,27 @@ function _pardiso_sequential_LODF!(
     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},
@@ -193,7 +214,11 @@ function _calculate_LODF_matrix_MKLPardiso(
     end
     lodf_t = zeros(linecount, linecount)
     A = SparseArrays.sparse(m_I, m_I, m_V)
-    _pardiso_sequential_LODF!(lodf_t, A, ptdf_denominator_t)
+    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

From c26a18c9dd831e93302800457a45a84271c7f98f Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 19:06:50 -0600
Subject: [PATCH 27/56] increase LODF chunk size

---
 src/definitions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/definitions.jl b/src/definitions.jl
index 1501c6ff..167eaf6a 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -6,4 +6,4 @@ const GiB = MiB * KiB
 const MAX_CACHE_SIZE_MiB = 100
 const ROW_PERSISTENT_CACHE_WARN = 1 * GiB
 
-DEFAULT_LODF_CHUNK_SIZE = 10_000
+DEFAULT_LODF_CHUNK_SIZE = 30_000

From 9a9baddbe5d1709d2bf7fd1b448be4ef5ba2628d Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 21:00:24 -0600
Subject: [PATCH 28/56] reduce chunk size

---
 src/definitions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/definitions.jl b/src/definitions.jl
index 167eaf6a..be3717f9 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -6,4 +6,4 @@ const GiB = MiB * KiB
 const MAX_CACHE_SIZE_MiB = 100
 const ROW_PERSISTENT_CACHE_WARN = 1 * GiB
 
-DEFAULT_LODF_CHUNK_SIZE = 30_000
+DEFAULT_LODF_CHUNK_SIZE = 20_000

From 3037b531bd17f2fc5c268a8dbc7470a92c347283 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 21:44:06 -0600
Subject: [PATCH 29/56] try using copyto

---
 src/ptdf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 558880b3..0d9beba8 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -291,13 +291,13 @@ function _calculate_PTDF_matrix_MKLPardiso(
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
-        PTDFm_t[valid_ix, :] = full_BA
+        copyto!(PTDFm_t[valid_ix, :], full_BA)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
-        PTDFm_t[valid_ix, :] = full_BA
+        copyto!(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)

From a5b6e4c438a118b09c921a2439b1db092275d96c Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 21:52:20 -0600
Subject: [PATCH 30/56] use dot

---
 src/ptdf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 0d9beba8..0e0260a2 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -291,13 +291,13 @@ function _calculate_PTDF_matrix_MKLPardiso(
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
-        copyto!(PTDFm_t[valid_ix, :], 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"
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
-        copyto!(PTDFm_t[valid_ix, :], 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)

From ee06bf84994f35c9f02a7516b230c27a9e0d26fa Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:05:33 -0600
Subject: [PATCH 31/56] not call release all

---
 src/ptdf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 0e0260a2..bc52d227 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -292,13 +292,13 @@ function _calculate_PTDF_matrix_MKLPardiso(
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         PTDFm_t[valid_ix, :] .= full_BA
-        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        #Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         PTDFm_t[valid_ix, :] .= full_BA
-        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        #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)

From 2e1eb70a50a3e44feed8358511aaea9a253c4df4 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:12:17 -0600
Subject: [PATCH 32/56] call pardiso again after release

---
 src/ptdf_calculations.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index bc52d227..cd1fd0b1 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -292,13 +292,15 @@ function _calculate_PTDF_matrix_MKLPardiso(
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         PTDFm_t[valid_ix, :] .= full_BA
-        #Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         PTDFm_t[valid_ix, :] .= full_BA
-        #Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)

From 643f7d3071b2f3b2b5f0ad9c647d59aca4bc4f52 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:18:53 -0600
Subject: [PATCH 33/56] one more try

---
 src/ptdf_calculations.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index cd1fd0b1..5c7ff52b 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -290,14 +290,14 @@ function _calculate_PTDF_matrix_MKLPardiso(
             "Distibuted slack is not supported for systems with multiple reference buses.",
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
-        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        Pardiso.pardiso(ps, similar(full_BA), ABA, full_BA)
         PTDFm_t[valid_ix, :] .= full_BA
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        Pardiso.pardiso(ps, Float64[], ABA, full_BA)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        Pardiso.pardiso(ps, similar(full_BA), ABA, full_BA)
         PTDFm_t[valid_ix, :] .= full_BA
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)

From 62370e4fda38a48e221d529ca85cc9fa4b1c612d Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:28:43 -0600
Subject: [PATCH 34/56] change iparm 12

---
 src/lodf_calculations.jl | 8 ++++++++
 src/ptdf_calculations.jl | 2 ++
 2 files changed, 10 insertions(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index cf849d0b..18720dc6 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -183,6 +183,14 @@ function _pardiso_single_LODF!(
     linecount = size(lodf_t, 1)
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.pardisoinit(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)
+    Pardiso.set_iparm!(ps, 12, 1)
     #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     Pardiso.pardiso(
         ps,
diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 5c7ff52b..0208be9e 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -270,6 +270,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
     ABA = calculate_ABA_matrix(A, BA, ref_bus_positions)
     # Here add the subnetwork detection
     ps = Pardiso.MKLPardisoSolver()
+    Pardiso.pardisoinit(ps)
     # Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     defaults = Pardiso.get_iparms(ps)
     Pardiso.set_iparm!(ps, 1, 1)
@@ -279,6 +280,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
     Pardiso.set_iparm!(ps, 2, 2)
     Pardiso.set_iparm!(ps, 59, 2)
     Pardiso.set_iparm!(ps, 6, 1)
+    Pardiso.set_iparm!(ps, 12, 1)
 
     # inizialize matrices for evaluation
     valid_ix = setdiff(1:buscount, ref_bus_positions)

From 8e3708fca8ddace1a859072f1306700d70025988 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:30:30 -0600
Subject: [PATCH 35/56] call pardiso after release

---
 src/lodf_calculations.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 18720dc6..1f024918 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -171,6 +171,11 @@ function _pardiso_sequential_LODF!(
         i_count = edge + 1
     end
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+    Pardiso.pardiso(
+        ps,
+        A,
+        Float64[],
+    )
     return
 end
 
@@ -199,6 +204,11 @@ function _pardiso_single_LODF!(
         ptdf_denominator_t,
     )
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+    Pardiso.pardiso(
+        ps,
+        A,
+        Float64[],
+    )
     return
 end
 

From 435131533ce22f905b8b6f639ba28ee28ff6644f Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:41:34 -0600
Subject: [PATCH 36/56] try with deepcopy

---
 src/ptdf_calculations.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 0208be9e..948c1af0 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -292,15 +292,15 @@ function _calculate_PTDF_matrix_MKLPardiso(
             "Distibuted slack is not supported for systems with multiple reference buses.",
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
-        Pardiso.pardiso(ps, similar(full_BA), ABA, full_BA)
-        PTDFm_t[valid_ix, :] .= full_BA
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        PTDFm_t[valid_ix, :] .= deepcopy(full_BA)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, Float64[], ABA, full_BA)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        Pardiso.pardiso(ps, similar(full_BA), ABA, full_BA)
-        PTDFm_t[valid_ix, :] .= full_BA
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        PTDFm_t[valid_ix, :] .= deepcopy(full_BA)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
         slack_array = dist_slack / sum(dist_slack)

From a14f03245b04a5f5e0689b3c089c962e46c29632 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:56:39 -0600
Subject: [PATCH 37/56] try to multiply

---
 src/ptdf_calculations.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 948c1af0..a043d0ba 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -293,6 +293,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        full_BA*1.0
         PTDFm_t[valid_ix, :] .= deepcopy(full_BA)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, Float64[], ABA, full_BA)

From e56491b77541da89aba4fbd165f0a865c5379ecf Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 22:57:16 -0600
Subject: [PATCH 38/56] change order of operations

---
 src/ptdf_calculations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index a043d0ba..476fbe45 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -293,10 +293,10 @@ function _calculate_PTDF_matrix_MKLPardiso(
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        Pardiso.pardiso(ps)
         full_BA*1.0
         PTDFm_t[valid_ix, :] .= deepcopy(full_BA)
-        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-        Pardiso.pardiso(ps, Float64[], ABA, full_BA)
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"

From 77b36f103deeda9b4fc5d2d0d4559afe26fc9e8a Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 23:06:11 -0600
Subject: [PATCH 39/56] dont store in full_BA

---
 src/ptdf_calculations.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 476fbe45..6d1d9618 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -279,7 +279,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
     end
     Pardiso.set_iparm!(ps, 2, 2)
     Pardiso.set_iparm!(ps, 59, 2)
-    Pardiso.set_iparm!(ps, 6, 1)
+    #Pardiso.set_iparm!(ps, 6, 1)
     Pardiso.set_iparm!(ps, 12, 1)
 
     # inizialize matrices for evaluation
@@ -292,11 +292,11 @@ function _calculate_PTDF_matrix_MKLPardiso(
             "Distibuted slack is not supported for systems with multiple reference buses.",
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
-        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        tmp = similar(full_BA)
+        Pardiso.pardiso(ps, tmp, ABA, full_BA)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps)
-        full_BA*1.0
-        PTDFm_t[valid_ix, :] .= deepcopy(full_BA)
+        PTDFm_t[valid_ix, :] = tmp
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"

From fd84136a23ebd14f04dfe9a5ac2723a31c82de77 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 23:13:44 -0600
Subject: [PATCH 40/56] add logging

---
 src/ptdf_calculations.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 6d1d9618..0894b8d5 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -293,9 +293,13 @@ function _calculate_PTDF_matrix_MKLPardiso(
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
         tmp = similar(full_BA)
+        @error "solve call"
         Pardiso.pardiso(ps, tmp, ABA, full_BA)
+        @error "solve done"
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+        @error "release call"
         Pardiso.pardiso(ps)
+        @error "release done"
         PTDFm_t[valid_ix, :] = tmp
         return PTDFm_t
     elseif length(dist_slack) == buscount

From 90cf3242e09656c9aebce463feb330c9d4a51249 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 23:18:34 -0600
Subject: [PATCH 41/56] fix to LODF

---
 src/lodf_calculations.jl | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 1f024918..1e953a55 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -189,12 +189,12 @@ function _pardiso_single_LODF!(
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.pardisoinit(ps)
     Pardiso.set_iparm!(ps, 1, 1)
+    defaults = Pardiso.get_iparms(ps)
     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)
     Pardiso.set_iparm!(ps, 12, 1)
     #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     Pardiso.pardiso(
@@ -204,11 +204,7 @@ function _pardiso_single_LODF!(
         ptdf_denominator_t,
     )
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-    Pardiso.pardiso(
-        ps,
-        A,
-        Float64[],
-    )
+    Pardiso.pardiso(ps)
     return
 end
 

From 53b09f8bf7415404fbf186502c41a22429991f20 Mon Sep 17 00:00:00 2001
From: jd-lara <jdlara@berkeley.edu>
Date: Wed, 27 Sep 2023 23:20:13 -0600
Subject: [PATCH 42/56] revert changes

---
 src/ptdf_calculations.jl | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 0894b8d5..7db51cd7 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -279,7 +279,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
     end
     Pardiso.set_iparm!(ps, 2, 2)
     Pardiso.set_iparm!(ps, 59, 2)
-    #Pardiso.set_iparm!(ps, 6, 1)
+    Pardiso.set_iparm!(ps, 6, 1)
     Pardiso.set_iparm!(ps, 12, 1)
 
     # inizialize matrices for evaluation
@@ -292,22 +292,17 @@ function _calculate_PTDF_matrix_MKLPardiso(
             "Distibuted slack is not supported for systems with multiple reference buses.",
         )
     elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
-        tmp = similar(full_BA)
-        @error "solve call"
-        Pardiso.pardiso(ps, tmp, ABA, full_BA)
-        @error "solve done"
+        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        PTDFm_t[valid_ix, :] = full_BA
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-        @error "release call"
         Pardiso.pardiso(ps)
-        @error "release done"
-        PTDFm_t[valid_ix, :] = tmp
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
         Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
-        PTDFm_t[valid_ix, :] .= deepcopy(full_BA)
+        PTDFm_t[valid_ix, :] = full_BA
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-        Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
+        Pardiso.pardiso(ps)
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)

From b240879d724255a109249dc9e26396d8d13602ae Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 01:38:20 -0600
Subject: [PATCH 43/56] update pariso params

---
 src/lodf_calculations.jl | 17 ++++++-----------
 src/ptdf_calculations.jl |  8 ++++++--
 2 files changed, 12 insertions(+), 13 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 1e953a55..10e98dc0 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -135,10 +135,10 @@ function _pardiso_sequential_LODF!(
     chunk_size::Int = DEFAULT_LODF_CHUNK_SIZE,
 )
     linecount = size(lodf_t, 1)
+    @assert LinearAlgebra.ishermitian(A)
     ps = Pardiso.MKLPardisoSolver()
-    Pardiso.pardisoinit(ps)
+    Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
     #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
-
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
     Pardiso.pardiso(
         ps,
@@ -157,7 +157,7 @@ function _pardiso_sequential_LODF!(
     i_count = 1
     tmp = zeros(Float64, linecount, chunk_size)
     while i_count <= linecount
-        edge = min(i_count + chunk_size -  1, linecount)
+        edge = min(i_count + chunk_size - 1, linecount)
         if linecount - edge <= 0
             tmp = tmp[:, 1:(edge - i_count + 1)]
         end
@@ -171,11 +171,7 @@ function _pardiso_sequential_LODF!(
         i_count = edge + 1
     end
     Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-    Pardiso.pardiso(
-        ps,
-        A,
-        Float64[],
-    )
+    Pardiso.pardiso(ps)
     return
 end
 
@@ -183,10 +179,10 @@ 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)
+    @assert LinearAlgebra.ishermitian(A)
     ps = Pardiso.MKLPardisoSolver()
+    Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
     Pardiso.pardisoinit(ps)
     Pardiso.set_iparm!(ps, 1, 1)
     defaults = Pardiso.get_iparms(ps)
@@ -208,7 +204,6 @@ function _pardiso_single_LODF!(
     return
 end
 
-
 function _calculate_LODF_matrix_MKLPardiso(
     a::SparseArrays.SparseMatrixCSC{Int8, Int},
     ptdf::Matrix{Float64},
diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 7db51cd7..ae85f017 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -268,10 +268,11 @@ function _calculate_PTDF_matrix_MKLPardiso(
     buscount = size(BA, 1)
 
     ABA = calculate_ABA_matrix(A, BA, ref_bus_positions)
-    # Here add the subnetwork detection
+    @assert LinearAlgebra.issymmetric(ABA)
     ps = Pardiso.MKLPardisoSolver()
+    Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
     Pardiso.pardisoinit(ps)
-    # Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    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])
@@ -281,6 +282,9 @@ function _calculate_PTDF_matrix_MKLPardiso(
     Pardiso.set_iparm!(ps, 59, 2)
     Pardiso.set_iparm!(ps, 6, 1)
     Pardiso.set_iparm!(ps, 12, 1)
+    Pardiso.set_iparm!(ps, 11, 0)
+    Pardiso.set_iparm!(ps, 13, 0)
+    Pardiso.set_iparm!(ps, 32, 1)
 
     # inizialize matrices for evaluation
     valid_ix = setdiff(1:buscount, ref_bus_positions)

From 03acf047ed61db28281ff133da3fe0dbd0e41289 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 01:47:12 -0600
Subject: [PATCH 44/56] remove messaging

---
 src/ptdf_calculations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index ae85f017..a2479a40 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -272,7 +272,7 @@ function _calculate_PTDF_matrix_MKLPardiso(
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
     Pardiso.pardisoinit(ps)
-    Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    # 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])

From 6fa823c64298cd15e33456a2ccb8de9878c2e8bb Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 02:54:52 -0600
Subject: [PATCH 45/56] use mkl linear algebra

---
 Project.toml                | 2 ++
 src/PowerNetworkMatrices.jl | 5 +++--
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 768294c6..4bf4c6bc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -17,6 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DocStringExtensions = "~0.8, ~0.9"
 HDF5 = "0.17"
 InfrastructureSystems = "^1.20"
+MKL = "0.6"
 KLU = "~0.4"
 Pardiso = "~0.5.4"
 PowerSystems = "3"
diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl
index abb3210d..53112ed4 100644
--- a/src/PowerNetworkMatrices.jl
+++ b/src/PowerNetworkMatrices.jl
@@ -19,6 +19,7 @@ export VirtualPTDF
 export Ybus
 
 using DocStringExtensions
+import MKL
 import InfrastructureSystems
 import PowerSystems
 import PowerSystems: ACBusTypes
@@ -31,9 +32,9 @@ import SparseArrays: rowvals, nzrange
 import HDF5
 import KLU: klu
 import KLU
-import LinearAlgebra: LAPACK.getri!, LAPACK.getrf!, BLAS.gemm
-import LinearAlgebra: ldiv!, mul!, I, dot
 import LinearAlgebra
+import LinearAlgebra: BLAS.gemm
+import LinearAlgebra: ldiv!, mul!, I, dot
 import Pardiso
 
 @template DEFAULT = """

From 10781b8a56a80a499e23f54663968d954dac3032 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 11:52:35 -0600
Subject: [PATCH 46/56] update pardiso calls

---
 src/lodf_calculations.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 10e98dc0..cc44d64a 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -134,10 +134,11 @@ function _pardiso_sequential_LODF!(
     ptdf_denominator_t::Matrix{Float64},
     chunk_size::Int = DEFAULT_LODF_CHUNK_SIZE,
 )
+    @info "Line Count too large for single compute using Pardiso. Employing Sequential Calulations using a chunk_size=$(chunk_size)"
     linecount = size(lodf_t, 1)
     @assert LinearAlgebra.ishermitian(A)
     ps = Pardiso.MKLPardisoSolver()
-    Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
+    Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
     #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
     Pardiso.pardiso(
@@ -170,8 +171,6 @@ function _pardiso_sequential_LODF!(
         lodf_t[:, i_count:edge] .= tmp
         i_count = edge + 1
     end
-    Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
-    Pardiso.pardiso(ps)
     return
 end
 

From 54a644c7a33fb9787da2fbf9102cd810b1e8c321 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 12:46:35 -0600
Subject: [PATCH 47/56] change parameters pardiso LODF

---
 src/lodf_calculations.jl | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index cc44d64a..575d69ed 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -139,6 +139,19 @@ function _pardiso_sequential_LODF!(
     @assert LinearAlgebra.ishermitian(A)
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
+    Pardiso.pardisoinit(ps)
+    # 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_iparm!(ps, 12, 1)
+    Pardiso.set_iparm!(ps, 11, 0)
+    Pardiso.set_iparm!(ps, 13, 0)
+    Pardiso.set_iparm!(ps, 32, 1)
     #Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
     Pardiso.pardiso(

From 5780ac9fab0abfa45ea441d3e1afe5fffdff7536 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 13:15:39 -0600
Subject: [PATCH 48/56] readd release

---
 src/lodf_calculations.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 575d69ed..ff58a44e 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -184,6 +184,8 @@ function _pardiso_sequential_LODF!(
         lodf_t[:, i_count:edge] .= tmp
         i_count = edge + 1
     end
+    Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
+    Pardiso.pardiso(ps)
     return
 end
 

From a471447339a56d0ed3b0a841a7cda4748a8fe4b1 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 15:05:24 -0600
Subject: [PATCH 49/56] reduce chunk size

---
 src/definitions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/definitions.jl b/src/definitions.jl
index be3717f9..1501c6ff 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -6,4 +6,4 @@ const GiB = MiB * KiB
 const MAX_CACHE_SIZE_MiB = 100
 const ROW_PERSISTENT_CACHE_WARN = 1 * GiB
 
-DEFAULT_LODF_CHUNK_SIZE = 20_000
+DEFAULT_LODF_CHUNK_SIZE = 10_000

From 96975cce42c1368b1c9263d07a865a5fb154e853 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 15:59:14 -0600
Subject: [PATCH 50/56] bump chunk size

---
 src/definitions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/definitions.jl b/src/definitions.jl
index 1501c6ff..819aaf7a 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -6,4 +6,4 @@ const GiB = MiB * KiB
 const MAX_CACHE_SIZE_MiB = 100
 const ROW_PERSISTENT_CACHE_WARN = 1 * GiB
 
-DEFAULT_LODF_CHUNK_SIZE = 10_000
+DEFAULT_LODF_CHUNK_SIZE = 15_000

From d2ed1ac0a9443a34bdc8b68a10ea5a105f966928 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 16:22:27 -0600
Subject: [PATCH 51/56] bump chunk size

---
 src/definitions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/definitions.jl b/src/definitions.jl
index 819aaf7a..b4d071a0 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -6,4 +6,4 @@ const GiB = MiB * KiB
 const MAX_CACHE_SIZE_MiB = 100
 const ROW_PERSISTENT_CACHE_WARN = 1 * GiB
 
-DEFAULT_LODF_CHUNK_SIZE = 15_000
+DEFAULT_LODF_CHUNK_SIZE = 18_000

From 08f168dadf215a0714e081446350a2ead3c1210d Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 28 Sep 2023 19:59:45 -0600
Subject: [PATCH 52/56] use openblas again

---
 Project.toml                | 4 ++--
 src/PowerNetworkMatrices.jl | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4bf4c6bc..b1b90efa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+# MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -18,7 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DocStringExtensions = "~0.8, ~0.9"
 HDF5 = "0.17"
 InfrastructureSystems = "^1.20"
-MKL = "0.6"
+# MKL = "0.6"
 KLU = "~0.4"
 Pardiso = "~0.5.4"
 PowerSystems = "3"
diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl
index 53112ed4..a90ec813 100644
--- a/src/PowerNetworkMatrices.jl
+++ b/src/PowerNetworkMatrices.jl
@@ -19,7 +19,7 @@ export VirtualPTDF
 export Ybus
 
 using DocStringExtensions
-import MKL
+#import MKL
 import InfrastructureSystems
 import PowerSystems
 import PowerSystems: ACBusTypes

From 8fb402b5e6050b6cbffaf6a601dd60ac098090c2 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Fri, 29 Sep 2023 08:25:14 -0600
Subject: [PATCH 53/56] se lapack directly

---
 src/PowerNetworkMatrices.jl |  3 ++-
 src/ptdf_calculations.jl    | 27 +++++++++++++++++++++------
 2 files changed, 23 insertions(+), 7 deletions(-)

diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl
index a90ec813..1546cd00 100644
--- a/src/PowerNetworkMatrices.jl
+++ b/src/PowerNetworkMatrices.jl
@@ -19,7 +19,7 @@ export VirtualPTDF
 export Ybus
 
 using DocStringExtensions
-#import MKL
+import MKL
 import InfrastructureSystems
 import PowerSystems
 import PowerSystems: ACBusTypes
@@ -35,6 +35,7 @@ import KLU
 import LinearAlgebra
 import LinearAlgebra: BLAS.gemm
 import LinearAlgebra: ldiv!, mul!, I, dot
+import LinearAlgebra: LAPACK.getrf!, LAPACK.getrs!
 import Pardiso
 
 @template DEFAULT = """
diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index a2479a40..03e62a08 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -123,19 +123,17 @@ function _calculate_PTDF_matrix_KLU(
     # inizialize matrices for evaluation
     valid_ix = setdiff(1:buscount, ref_bus_positions)
     PTDFm_t = zeros(buscount, linecount)
-
+    copyto!(PTDFm_t, BA)
     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
-        copyto!(PTDFm_t, BA)
         PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :])
         PTDFm_t[collect(ref_bus_positions), :] .= 0.0
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        copyto!(PTDFm_t, BA)
         PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :])
         PTDFm_t[collect(ref_bus_positions), :] .= 0.0
         slack_array = dist_slack / sum(dist_slack)
@@ -172,6 +170,19 @@ 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.
 
@@ -198,17 +209,21 @@ function _calculate_PTDF_matrix_DENSE(
     valid_ixs = setdiff(1:buscount, ref_bus_positions)
     ABA = Matrix(calculate_ABA_matrix(A, BA, ref_bus_positions))
     PTDFm_t = zeros(buscount, linecount)
-
+    (ABA, bipiv, binfo) = getrf!(ABA)
+    _binfo_check(binfo)
+    BA = Matrix(BA[valid_ixs, :])
     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, :] = ABA \ BA[valid_ixs, :]
+        getrs!('N', ABA, bipiv, BA)
+        PTDFm_t[valid_ixs, :] = BA
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        PTDFm_t[valid_ixs, :] = ABA \ BA[valid_ixs, :]
+        getrs!('N', ABA, bipiv, BA)
+        PTDFm_t[valid_ixs, :] = BA
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t -

From 451381c614778fe0a5f43d76349a1b4844a43fab Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Fri, 29 Sep 2023 08:25:20 -0600
Subject: [PATCH 54/56] add MKL again

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index b1b90efa..4bf4c6bc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-# MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -18,7 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DocStringExtensions = "~0.8, ~0.9"
 HDF5 = "0.17"
 InfrastructureSystems = "^1.20"
-# MKL = "0.6"
+MKL = "0.6"
 KLU = "~0.4"
 Pardiso = "~0.5.4"
 PowerSystems = "3"

From 9eea1369abb46a9aea11a58742a77c6e84b8dbe2 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Fri, 29 Sep 2023 08:56:50 -0600
Subject: [PATCH 55/56] update lodf LAPACK code

---
 src/lodf_calculations.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index ff58a44e..baaf79ff 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -123,9 +123,11 @@ function _calculate_LODF_matrix_DENSE(
             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
+    (mV, bipiv, binfo) = getrf!(Matrix(LinearAlgebra.diagm(m_V)))
+    _binfo_check(binfo)
+    getrs!('N', mV, bipiv, ptdf_denominator_t)
+    ptdf_denominator_t[LinearAlgebra.diagind(ptdf_denominator_t)] .= -1.0
+    return ptdf_denominator_t
 end
 
 function _pardiso_sequential_LODF!(

From 74e093fb8003dd30190ed66af7720ff980fde0e2 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Fri, 29 Sep 2023 09:33:59 -0600
Subject: [PATCH 56/56] update lodf klu

---
 src/lodf_calculations.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index baaf79ff..98f12076 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -80,9 +80,9 @@ function _calculate_LODF_matrix_KLU(
         end
     end
     Dem_LU = klu(SparseArrays.sparse(m_I, m_I, m_V))
-    lodf_t = Dem_LU \ ptdf_denominator
-    lodf_t[SparseArrays.diagind(lodf_t)] .= -1.0
-    return lodf_t
+    KLU.solve!(Dem_LU, ptdf_denominator)
+    ptdf_denominator[SparseArrays.diagind(ptdf_denominator)] .= -1.0
+    return ptdf_denominator
 end
 
 function _calculate_LODF_matrix_KLU(